Benthic communities under anthropogenic pressure show resilience across the Quaternary

The Southeast Pacific is characterized by rich upwelling systems that have sustained and been impacted by human groups for at least 12 ka. Recent fishing and aquaculture practices have put a strain on productive coastal ecosystems from Tongoy Bay, in north-central Chile. We use a temporal baseline to determine whether potential changes to community structure and composition over time are due to anthropogenic factors, natural climatic variations or both. We compiled a database (n = 33 194) with mollusc species abundances from the Mid-Pleistocene, Late Pleistocene, Holocene, dead shell assemblages and live-sampled communities. Species richness was not significantly different, neither were diversity and evenness indices nor rank abundance distributions. There is, however, an increase in relative abundance for the cultured scallop Argopecten, while the previously dominant clam Mulinia is locally very rare. Results suggest that impacts from both natural and anthropogenic stressors need to be better understood if benthic resources are to be preserved. These findings provide the first Pleistocene temporal baseline for the south Pacific that shows that this highly productive system has had the ability to recover from past alterations, suggesting that if monitoring and management practices continue to be implemented, moderately exploited communities from today have hopes for recovery.

The Southeast Pacific is characterized by rich upwelling systems that have sustained and been impacted by human groups for at least 12 ka. Recent fishing and aquaculture practices have put a strain on productive coastal ecosystems from Tongoy Bay, in north-central Chile. We use a temporal baseline to determine whether potential changes to community structure and composition over time are due to anthropogenic factors, natural climatic variations or both. We compiled a database (n = 33 194) with mollusc species abundances from the Mid-Pleistocene, Late Pleistocene, Holocene, dead shell assemblages and live-sampled communities. Species richness was not significantly different, neither were diversity and evenness indices nor rank abundance distributions. There is, however, an increase in relative abundance for the cultured scallop Argopecten, while the previously dominant clam Mulinia is locally very rare. Results suggest that impacts from both natural and anthropogenic stressors need to be better understood if benthic resources are to be preserved. These findings provide the first Pleistocene temporal baseline for the south Pacific that shows that this highly productive system has had the ability to recover from past alterations, suggesting that if monitoring and management practices continue to be implemented, moderately exploited communities from today have hopes for recovery.

Introduction
Humans have been exploiting and modifying marine environments around the world for ca 160 ka [1][2][3][4][5][6]. This pressure has increased over the last centuries and decades, leading to fundamental changes in the structure and functioning of these environments [3,4,7,8]. Classical studies incorporating historical, archaeological and palaeontological data have shown that overfishing, in particular, has driven once-abundant large pelagic predators to extinction [5] and has also dramatically altered coastal ecosystems [5,7]. Even if different human groups have been making use of coastal resources for millennia [1,2,[9][10][11], evidence shows that prehistoric human populations did not exert such a strong negative pressure, in comparison to the one marine systems have been suffering in the most recent industrial times [7,10,12,13]. In addition to overfishing, coastal areas are subject to habitat modification, deviation of watercourses, runoff of pollutants, aquaculture, and among others [7,12,14,15]. For example, the damming of the Colorado river caused the collapse of once very abundant populations of the clam Mulinia coloradoensis, leading to a dramatic drop in productivity [15], changes to trophic structure [16] and a reduction in carbon emission in the river basin [17]. Ever-growing demands for more food have led to poor aquaculture practices that bring about changes in community structure, ecosystem function, eutrophication and outbreaks of disease [5,18]. Yet, despite considerable efforts to summarize global trends and patterns in exploitation and degradation of coastal areas, most of the available information is limited to the Northern Hemisphere (but see [19] for Tasmania). Our present paradigm is thus lacking information from other coastal environments, also under anthropogenic strain, that can provide insights for a more holistic understanding and potentially a more holistic approach towards remediation and conservation. In this context, coastal marine environments from the Southeast Pacific, where there are highly productive fisheries related to the Humboldt Current System [13,20] are key pieces to add to the puzzle.
The Southeast Pacific is characterized by important upwelling systems along the Chilean and Peruvian coasts [20][21][22], that have sustained and been impacted by human groups for at least 12 ka [11,[23][24][25][26][27]. One of the most productive of these upwellings is located in the north-central region of Chile, near Tongoy (30°12 S-71°34 W, [20,21,28], figure 1). Owing to its high productivity this bay has seen the development of small-scale benthic fisheries and of an aquaculture regime of the scallop Argopecten purpuratus [29,30]. The uncontrolled exploitation of the scallop since 1945 led to the collapse of this fishery [25,29,30] but, given the economic importance of this resource, an aquaculture regime was implemented in 1988 [29,30]. The area destined for scallop culture is 54% of the surface of Tongoy Bay (1900 out of 3500 ha), and of the remaining 1500 ha, 100 ha are managed by local fishermen as a 'Management Exploitation Area for Benthic Resources' (AMERB) [30].
The modern benthic assemblages from Tongoy Bay have been extensively sampled and studied from a trophic network approach [31][32][33][34]. These contributions show that the uncontrolled fishing and aquaculture regimes have not only directly impacted the exploited species, but also the structure and functioning of the ecosystem they are a part of [31]. These trophic network studies are incredibly valuable to understand changes to the system in the last 30 years; however, information prior to the 1980s is lacking [31]. In order to assess the magnitude of these changes, it is necessary to use the temporal perspective provided by a Conservation Palaeobiology approach, to compare the current community to a pre-human-impact state [7,15,19,[35][36][37][38][39][40][41][42]. Tongoy Bay has dated marine terraces from the Pleistocene and Holocene (HOL) [43][44][45][46] Here, we seek to use a temporal baseline to contextualize the effect of artisanal fishing and a recent aquaculture regime to a productive coastal ecosystem from north-central Chile. We use palaeontological data from the Pleistocene and Holocene terraces together with recent ecological sampling to quantify community structure and composition over time, to determine whether there are changes and whether these are due to human-related pressures, natural climatic variations or both. We hypothesize that anthropogenic changes have altered community structure and composition over time, and that these will be significantly different in previous 'states' of the community. The results from this study will provide one of the first temporal baselines of this kind for the Southern Hemisphere, and in particular for a highly productive upwelling system from the Southeast Pacific.

Compilation of datasets
We compiled a database using published literature and collected fossil samples from Tongoy Bay. For the Mid-Pleistocene, 22 samples were used (all collected); for the Late Pleistocene, 26 (all collected); for the Holocene, 40 (all collected); for Dead Assemblages, 11 (three collected and eight from the published literature) and for the Live samples, we used a live-collected dataset from [47] and a 2012 live-collected dataset (partly published in [31], table 1). All the fossil samples correspond to bulk-collected unlithified samples comparable to the unlithified Dead Assemblages and the Live samples. Quaternary terraces of Tongoy have been dated by previous studies, showing a strong correlation between height above sea level and age [45,46]; hence, fossil sites were assigned to different interglacial stages from their height above sea level. Species from each time bin were pooled and labelled as follows: Mid-Pleistocene (samples from MIS7 and MIS9 where pooled together and labelled MP), Late Pleistocene, Holocene, Dead Assemblage (DA), Live community (LIVE). For some analyses, we pooled together MP and LP into 'Pleistocene' and DA and LIVE into 'Modern' (see details below).
For the fossil samples, shells from each time bin were identified to species level using published literature [48][49][50][51]. After species were identified, bivalve and gastropod individuals were counted. For bivalves, the number of individuals was calculated diving the total number of valves by two. The livecollected data from [47] (table 1) was collected using transects that were laid perpendicular to the bay. Samples were collected by divers using a quadrat at depths between 7 and 25 m. For the second, 2012 study, transects were also laid perpendicular to the bay, and quadrats were used to sample benthic organisms at 4, 8, 12 and 20 m depth. Only the data for bivalves and gastropods were used in order to have consistency with what is preserved in the fossil samples. In addition, the pooled live-collected data and the fossil data were checked for synonymy and updated taxonomy using the WoRMS taxon match online tool (http://www.marinespecies.org/). After the taxonomy was homogenized, we proceeded to separate species into categories depending if they were commercially exploited or non-exploited using the species inventories of the Servicio Nacional de Pesca [52]. The datasets supporting this article have been uploaded as part of the electronic supplementary material and are available on Dryad.

Diversity and abundance metrics
Species relative abundances per time bin (Pleistocene, Holocene and Modern) were calculated to determine whether the older time bin was a good predictor of the one that followed. That is, Pleistocene relative abundance (dependent variable) was regressed on Holocene relative abundance (independent variable), and Holocene relative abundance was regressed on Modern relative abundance. High agreement is indicated by species plotting along a 1 : 1 line in a bivariate plot of relative abundance [53,54]. To test whether the slopes of those regressions were significantly different from 1, we calculated the upper and lower confidence intervals. If a slope of 1 fell within the confidence intervals, we assumed that differences were not significant. The residuals for the regression with the total relative abundances were inspected to identify potential outliers. Species were considered outliers if they had a very high abundance in the dead assemblage and were absent or had low abundances in the living assemblage. For this analysis, we standardized Pleistocene, Holocene and Modern sample numbers by doing a rarefied subsample to the smallest sample number (n = 5176, for MP and LP pooled together). The rarefaction was done without replacement, using the 'rrarefy' function in the 'vegan' package [55], in the statistical programming language R [56].
Diversity indices were calculated for each time bin (MP, LP, HOL, DA and LIVE). For these analyses, we also standardized sample numbers by doing a rarefied subsample (without replacement) to the

Similarity metrics for the whole assemblage, exploited and non-exploited species
To test for similarity in species composition between the different time bins (MP, LP, HOL, DA and LIVE) Chao's Jaccard similarity index was used [59]. Chao's Jaccard index includes the effect of species that are shared but unseen (either because they are rare or because the samples that are being compared have substantial differences in size such as these live-dead assemblages). By accounting for unseen species, this estimator is less biased than the classic Jaccard index that is only based on the presenceabsence data [59]. The Spearman rank-order correlation of species relative abundance was also used as an indicator of similarity between the different time bins [60]. Chao's Jaccard similarity index and Spearman's rank-order correlation are typically plotted on bivariate plots to represent compositional and abundance similarity between assemblages. In this plot, sites located in the upper right-hand quadrant have the highest agreement and sites in the lower left-hand quadrant have the lowest agreement [53]. Indices were calculated with the 'diversity' and 'chao.jaccard' functions, in the 'vegan' [55] and 'fossil' [61] packages in the statistical programming language R [56].

Rank abundance distributions
Species rank abundance plots are also good descriptors of communities [62]. Several theories and models have been proposed to explain the different shape of rank abundance plots in communities (see [62] for a review). Here, we fit three of these models (Geometric series, Broken stick and Zipf) to the rank abundance orders of assemblages from different time bins (Pleistocene, Holocene and Modern) to determine the best-fit model for each dataset. The best model was chosen based on at least a two-point difference in Akaike Information Criterion (AIC). For this analysis, we also standardized Pleistocene, Holocene and Modern sample numbers by doing a rarefied subsample without replacement to n = 5176. We carried out these analyses with the 'rrarefy' and 'fitrad' functions in the 'vegan' [55] and 'sads' [63] packages in the statistical programming language R [56].

Abundance and body size of Argopecten and Mulinia through time
The relative abundance (i.e. proportion of the total individuals) of the scallop Argopecten purpuratus (Lamarck 1819) and the clam Mulinia edulis (King & Broderip 1831) were quantified per time bin. These two species were selected as they were very abundant in fossil samples, and are subject to fishing pressure and aquaculture (only Argopecten) nowadays. Changes in body size can be indicative of subsistence harvesting [64] and fishing pressure [65]. Species typically exhibit decreasing body size, indicative of an overexploitation of larger size classes [64,65]. The two most abundant species in Tongoy Bay are subject to different anthropogenic pressures as Argopecten is cultivated, but Mulinia is not. Thus, it is possible there are variations in size between the fossil and the dead samples for these species. To explore this, specimens from Argopecten (n = 135) and Mulinia (n = 3803) were measured to the nearest millimetre using a digital caliper, and their body size was calculated per time bin as the geometric mean of shell height and shell length [66].

Diversity and abundance metrics
The samples collected and compiled from the literature (n = 690, table 1) yielded a live assemblage with  18 161 individuals from 36 species, DAs with 1934 shells from 24 species, and fossil assemblages with  3242 individuals and 20 species from the HOL, 8189 and 24 from the LP, and 1668 and 19 from the MP  (table 1). The combined richness was of 62 species.
Univariate raw unstandardized diversity metrics are not different between time periods (table 1). Species richness is not significantly different (Kruskal-Wallis rank sum test, χ 2 = 4, p = 0.406), neither are Shannon's or Pielou's diversity and evenness indices (Kruskal-Wallis rank sum test, χ 2 = 4, p = 0.406 for both). When Pleistocene and Holocene samples are pooled together, linear models indicate that abundances of Pleistocene species are significant predictors of the Holocene species relative abundance (adjusted R 2 = 0.20, F = 5.83, p = 0.03, electronic supplementary material, figure S1a); however, pooled Holocene species are not significant predictors of Modern species relative abundances (adjusted R 2 = −0.05, F = 0.009, p = 0.92; electronic supplementary material, figure S1b), suggesting a shift in species composition and/or relative abundances in the live-collected samples. Confidence intervals indicate that slopes are slightly lower than 1 for the Pleistocene and Holocene (lower CI = 0.07, upper CI = 0.96), but the Holocene and Modern slope is significantly lower than 1 (lower CI = −0.48, upper CI = 0.53).

Similarity metrics for the whole assemblage, exploited and non-exploited species
Chao's Jaccard index for assemblage-level compositional similarity between time periods is higher than 0.64 for all comparisons (electronic supplementary material, table S1; figure 2a). This similarity between samples suggests stability in composition throughout the span of the Quaternary analysed. A visual inspection of the bivariate plot with Chao's Jaccard similarity index and Spearman's ρ (figure 2a) shows that the samples fall in the upper right-hand quadrant, indicating that live-dead agreement is high [53]. Significance for Spearman's rank correlations between samples is shown in the electronic supplementary material, table S1.
When species were subdivided into 'Exploited' and 'Non-exploited' categories it became evident that 'Exploited' species were driving the dissimilarity between live-collected and fossil samples (electronic supplementary material, table S1; figure 2b,c). The species classified as 'Exploited' (or 'exploitable' in the case of the fossil samples that are prior to documented human exploitation) make up a large proportion of the community. Moreover, 'Non-exploited' species tend to be less abundant, sometimes rare. Thus, it is not surprising that the data show a less clear pattern for the latter (figure 2b,c). Spearman's ρ is significant for only half of the comparisons in 'Exploited' and 'Non-exploited' (electronic supplementary material, table S1). A visual inspection of figure 2b,c shows that comparisons between fossils and the live-collected samples plot in the upper right-hand quadrant for 'Exploited', but differ for just 'Non-exploited'.

Rank abundance distributions
Species rank abundance distributions for Modern, Holocene and Pleistocene assemblages are best explained by the same model, as indicated by AIC (electronic supplementary material, table S2; figure 3a-c). The model with the strongest support is Zipf.

Abundance of Argopecten and Mulinia through time
The relative abundance of 'Exploited' and 'Non-exploited' species was calculated for each time bin (figure 4a,b). For the LIVE assemblage 21% of the individuals belong to exploited species and 50% of that corresponds to the scallop Argopecten. If these values are compared to those from the older fossil assemblages, the individuals that belong to exploited species make up between 40% (MP) and 92% (HOL) but Argopecten was responsible for less than a fifth of this, and the clam Mulinia for over 50% (figure 4b). Thus, locally harvested Argopecten shows an increase in relative abundance in LIVE compared with fossil assemblages, whereas Mulinia used to be very abundant in the recent past (over 80% of the 'exploitable' species) but has decreased in abundance and is rare in the LIVE assemblage ( figure 4b).     fishing and aquaculture pressure, in the last two decades the benthic community shows signs of resilience and recovery [31]. When results are considered for two of the most abundant species, however, this general positive outlook becomes questionable. There is an artificial increase in the relative abundance of Argopecten, while Mulinia is nowadays very rare, and only consistently found alive in neighbouring bays. The increase in Argopecten abundance is probably a reflection of aquaculture practices, yet we are unaware of any detrimental effects this practice could have had on Mulinia. The absence of this clam may be due to natural environmental changes in the region, local changes in freshwater input to the bay (see 'Reconstructing the regional history of Mulinia' in the Discussion), and/or to more recent anthropogenic impacts such as disease by parasites (see 'Recent local anthropogenic pressures on the benthic community' in the Discussion). Results indicate that natural and anthropogenic stressors probably impacted the populations at different times, with strong anthropogenic impacts being quite recent (less than 50 years). Findings from this study suggest that benthic communities from the south Pacific show strong resilience and recovery potential, but drivers of change and their impacts on key species need to be better understood if benthic resources are to be preserved in the near future. Moreover, even if it exceeds the scope of the manuscript, it is worth emphasizing that in order to reconstruct changes in the community over time a Conservation Palaeobiology approach is fundamental. This would require doing more detailed palaeoenvironmental reconstructions together with a refined stratigraphic framework.

Stability and recovery despite human pressure
The stability or resilience suggested by the mollusc data from Tongoy Bay is in agreement with global studies that also indicate high recovery potential in marine ecosystems [7,68]. For example, fossil marine faunas have shown stasis through time [69][70][71] while modern marine faunas seem to have high recovery potential, returning to conditions similar to the original one in decades [7,72,73]. More specifically for molluscs, it has been shown that the composition of bivalves and gastropods from reefs in the Bahamas had remained stable from the Pleistocene [74]. Other studies have nevertheless found remarkable differences between marine communities from the Pleistocene and Holocene with present-day assemblages [39,70]. For example, fossil molluscs from core samples collected in the Adriatic Sea had different composition, diversity and dominance than their recent counterparts, suggesting these communities remained relatively unchanged throughout glacial-interglacial cycles but shifted in composition in recent times probably due to anthropogenic impacts [39]. Similarly, coral communities from Barbados showed resilience and stability through the Pleistocene but dramatic shifts in composition and abundance were observed in recent coral species relative to fossil ones [70]. Results from these temporal studies can, however, be contingent on the scale of analyses [75]. Similar research on molluscs found that species composition was highly similar between Pleistocene and Modern samples when viewed at a regional scale, but similarity decreased when comparisons were done at a local scale [75]. In addition, preliminary results suggest strong differences in species composition between Modern and Pleistocene assemblages in other localities in north-central Chile [42].

Reconstructing the regional history of Mulinia
Despite the stability in composition and abundance in the mollusc community, and the lack of significant differences between the measured indices, Mulinia is locally very rare in Tongoy Bay and is only consistently found alive in neighbouring bays [67]. This marked change in the dominance of the clam can be a consequence of natural and anthropogenic causes. regional information, previous studies on congeneric species, and fisheries landings suggest that the disappearance of the clam was mainly due to natural factors. At a regional level, Mulinia was the most common bivalve in Pleistocene assemblages from northern Chile and Peru [44]. There are, however, notable differences in abundance and size for Mulinia between the MP and LP. This could be associated with the much higher uplift rates estimated for Tongoy Bay during the MP (in particular MIS 7 and MIS 9) compared with the LP [45]. These high uplift rates induced by tectonic processes could have led to a higher probability of coastal modification and habitat destruction, which could have impacted molluscan abundance and size. The species seems to have recovered in the LP. During the Holocene, however, the record of the Mulinia is very sparse, as the clam practically disappeared from northern Chile and is only found in one Holocene terrace in Peru (Michilla terrace around 7 ka [44]). Thus, it has been suggested that the dramatic shrinking of the species distribution occurred in the Holocene, possibly by a combination of palaeoceanographic circumstances together with biological phenomena [44]. More recent palaeoclimatic studies in northern Chile and Peru support this idea, given that aridity increased in the region during the Mid-Holocene [76][77][78][79]. Specific dates vary but most studies indicate that a displacement of the southern Westerlies (which control an N-S precipitation gradient) caused arid conditions between 7.7 and 4.2 ka [76][77][78][79]. This aridity could have indirectly affected Mulinia due to a decrease in freshwater input. Previous studies show that populations of the congeneric Mulinia coloradoensis from the Northern Hemisphere were decimated by a decrease in freshwater and nutrient input [15]. These changes caused by the damming of the Colorado river significantly decreased the productivity of the associated estuary and brought about detrimental ecosystem-level consequences [15][16][17]80]. Therefore, it is likely that Mulinia from Chile and Peru could have also thrived in estuary conditions during more humid weather. This idea is further supported by the presence of the species in estuaries in Los Lagos region, southern Chile.
Another local line of evidence suggesting that a freshwater input is important for the species comes from recent fisheries landings in the area (electronic supplementary material, figure S2). According to this information, the species was initially not exploited (first records are from 1994) due to low abundance but there was a boost in the population during 1997, a strong ENSO year [20], and a sharp decline in landings afterwards. The increase in rain that year brought water to previously dry riverbanks that discharge in Tongoy Bay, perhaps creating favourable conditions for the clam and leading to an increased catchment by the local fishermen towards the end of the year. Regardless of this isolated local event, if we consider together the evidence from regional climate changes and life-history information of congeneric species, we can suggest that the disappearance of the clam Mulinia was probably a regional phenomenon driven by natural changes in environmental conditions, rather than a consequence of negative human impacts.

Recent local anthropogenic pressures on the benthic community
Despite the strong evidence pointing to natural regional changes, human-induced aquaculture and fishing pressures might also be having a more recent detrimental effect. Studies looking at the status of the system in Tongoy Bay from a trophic network approach found that the overall health of the community improved from 1994 to 2012 [31]. This apparent recovery was brought about by release from fishing pressure due to the establishment of better management practices in the 1990s. Mulinia is, however, still not present in the living community, and the banks of the species in neighbouring bays are 97% infected with a trematode parasite that attacks soft tissue, the gills in particular [67]. No information is yet available on the parasite species or where it came from but its presence might be related to aquaculture practices in the bay [81]. For example, the scallop Argopecten has numerous described parasites and commensals [82] that could also use Mulinia as a host. Therefore, even if the major changes to species populations and abundance over time appear to have been natural, there are anthropogenic factors that presently pose a threat, might have not left a record yet and should thus continue to be closely monitored.

Caveats and closing remarks
So far, results indicate that the benthic environment is Tongoy Bay has been stable through the Quaternary and if changes occurred, the system showed resilience. Nevertheless, a caveat to consider is that this study was centred on the benthic mollusc community, limiting our understanding of large and/or pelagic organisms, for example. Previous research [3,5,83] has shown that large pelagic predators were decimated in tropical seas worldwide, leading to changes in functionality in the whole community. We are unaware if something like this happened to large pelagic fish in Tongoy Bay but the evidence at hand shows that mollusc communities are good surrogates for the benthic community [84], suggesting that our findings are at least a reliable representation of changes to this part of the community over time.
Here, we present what is to our knowledge the first Conservation Palaeobiology study for South America and for a highly productive upwelling system in the south Pacific Ocean. Our results suggest that the marine benthic community shows resilience and recovery from the Mid-Pleistocene onwards. An important decrease in relative abundance to the once dominant clam Mulinia seems to respond to natural climatic shifts in the Holocene. There are, however, indications that recent anthropogenic pressures may lead to unseen changes in the region as the clam is infested by parasites in neighbouring bays, and the scallop Argopecten has an artificially increased relative abundance probably due to aquaculture. These findings provide a temporal baseline that shows that this highly productive system has had the ability to recover from past alterations. Nonetheless, trophic network studies in Tongoy Bay [31] and in central Chile [85] suggest that fisheries can heavily modulate subtidal communities, with impacts extending even to non-harvested species [85]. Therefore, whether more recent human pressures will modify this long-time history of resilience remains to be seen. This temporal perspective provides an understanding of the variability displayed by this benthic marine system in the past, which is critical in order to manage for the future [36,37,86,87]. Studies of this kind are much needed to better comprehend recent changes to global marine communities where the literature is dominated by examples from the Northern Hemisphere and tropical environments. Assuming the system continues to behave as in the past, findings presented here suggest that if monitoring and management practices continue to be implemented, the moderately exploited communities from today have high hopes for recovery.
Data accessibility. Supporting data are available as supplementary online material (data, figures, tables and R code).