Seasonal patterns of biodiversity in Mediterranean coastal lagoons

Understanding and quantifying the seasonal patterns in biodiversity of phytobenthos, macro‐zoobenthos and fishes in Mediterranean coastal lagoons, and the species dependence upon environmental factors.

Shallow waters are mostly in close proximity to large human populations and are therefore affected by management-induced alterations to physical and chemical conditions (Airoldi & Beck, 2007) that may exacerbate global climate change effects on ecosystems (Anthony et al., 2009).
Easy to reach and rich in biodiversity, shallow waters not only provide a wide range of ecosystem services to citizens but they are also particularly suitable as open-air laboratories to study how biological diversity is affected by changing environmental conditions. Further, shallow waters can be data replete as they easily host large-scale studies measuring community biodiversity (Bintz, Nixon, Buckley, & Granger, 2003;Grenz et al., 2017) and a large number of environmental variables, allowing to investigate the driving processes of biodiversity variations. These systems are thus ideal to develop new avenues for predicting ecosystem responses at scales that are relevant to biodiversity management and conservation in the context of global climate change (Albouy et al., 2015;Petes, Howard, Helmuth, & Fly, 2014).
While our ability to provide future predictions is increasing at functional levels with the recent introduction of mechanistic individual trait-based bioenergetic models (e.g., Kearney & Porter, 2009;Sarà, Palmeri, Rinaldi, Montalto, & Helmuth, 2013;Sarà, Rinaldi, & Montalto, 2014;Sarà, Mangano, Johnson, & Mazzola, 2018), quantifying changes in community structure and biodiversity is more complex and challenging (Bellard et al., 2012;Best et al., 2015). On the one hand, machine-learning or deep-learning algorithms may provide a mathematical framework to model complex dynamics, although they may lack ecological realism. On the other hand, process-based models may provide exceptional realism, but the complexity required to model individual populations in complex communities may be overwhelming. Causal network techniques such as Bayesian graphical models, structural equation models or path analysis bridge the gap between the two approaches, providing ecological realism without the K E Y W O R D S artificial neural networks, biodiversity, climate change, community structure, confirmatory path analysis, fish, lagoon systems, phytobenthos, ridge regression, zoobenthos complexity of process-based models, but they are not particularly suited to model complex communities either. These models have been used, for example, to examine the relative weight of environmental factors and anthropogenic disturbances (Schoolmaster, Grace, Schweiger, Mitchell, & Guntenspergen, 2013), and the role of carnivore (Calcagno, Sun, Schmitz, & Loreau, 2011) or alien species on local community structure (McMahon, 2005). In marine systems, causal analysis has been used to model trophic web responses to sea surface temperature in the North Sea (Kirby & Beaugrand, 2009).
The "Stagnone di Marsala e Saline di Trapani e Paceco" is the largest coastal lagoon system in Sicily (Italy) and represents one of the central ecological corridors joining Africa and Europe, housing the largest stop-over sites for migrating avifauna in the central Mediterranean (Special Protection Area; EU Birds Directive 79/409/EEC). Here, we quantify biodiversity in different communities along the lagoon food web and clarify how their species respond to environmental factors.
Specifically, we aimed at evaluating the seasonal variations in biodiversity and at modelling the responses of individual species belonging to phytobenthic, zoobenthic and fish communities, to temperature, salinity, depth, organic and inorganic suspended matters. Finally, we exemplify a possible modelling approach, based on the relationships among communities and environmental variables, developing also future scenarios for the changes in coastal lagoon ecosystems driven by the environmental conditions expected under climate change.

| Study area
The study was carried out in Western Sicily (Southern Italy) in the "Stagnone di Marsala e Saline di Trapani e Paceco," one of the most important coastal lagoon systems of central Mediterranean Sea, encompassing 25 water bodies and representing an important hot spot of biodiversity for the Mediterranean Sea (Mannino, 2010;Mannino & Graziano, 2016;Mannino & Sarà, 2006;Mazzola et al., 2010;Vizzini, Sarà, Mateo, & Mazzola, 2003). The site, recog-  (Figure 1). The area, described in detail in Sarà, Leonardi, and Mazzola (1999), covers a total of 3,742 ha and includes both a large coastal pond system ("Saline di Trapani e Paceco") and the largest coastal lagoon in Southern Italy: the "Stagnone di Marsala" (Site of Community Importance since 1995, ITA010026, EU Habitat Directive 92/43/EEC). The entire area is subjected to a Mediterranean climate, with mild winters and the highest temperatures in summer, associated with the lowest precipitations.

| Field surveys and data collection
Sampling surveys were carried out in each season, during 2011-2012, in nine coastal ponds ( Figure 1 Pc and Zc communities were sampled manually over three quadrats of 20 cm 2 per pond in each season, whereas Fc by mean of three seine nets (12 × 1 m; 5 mm wings mesh, 1 mm coded mesh) deployed by hand, per pond in each season. Samples from both quadrats and nets were fixed in 70% ethanol and properly labelled. Once in the laboratory, samples were gently washed over sieves (0.5 mm mesh), and thereafter, Pc and Zc components were sorted. Pc and Zc were identified to the lower taxonomic level possible (Mannino & Sarà, 2006 and references therein; Mangano et al., 2014), and the total number of individuals per species was quantified. Environmental variables included in this study were measured near the bottom in triplicate per pond in each season. T, S and W were measured in the field near each quadrat with a multiprobe gauge (YellowSpring Ysi 556), whereas ISM and OSM were measured on water samples, collected nearby the quadrats using Niskin bottles, through drying and loss-on-ignition (450°C for 4 hr). Further details on the methods employed in environmental analyses can be found in Sarà et al. (1999).
Overall, the dataset comprised a number of 108 quadrats, seine nets and measures of each environmental variable (3 replicates × 9 basins × 4 seasons).

| Data analysis
The number of taxa and the relative abundances of Pc, Zc and Fc were measured and combined using biodiversity indices (Magurran, 2004): estimated species richness (S ACE ), Margalef (D Mg ), Shannon (H′), Brillouin (HB), Simpson (D) and Pielou (J′). Two indices for each biodiversity measure, that is, species richness (S ACE and D Mg ), diversity (H′ and HB) and dominance/evenness (D and J′) were employed.
The S ACE index was derived according to abundance-based coverage estimate (ACE) models. All the indices were calculated based on the absolute abundances, where applicable. Differences in biodiversity among seasons and communities were evaluated through a multivariate analysis of variance (MANOVA) using the biodiversity indices as dependent variables, the season and community variables as fixed factors and Pillai's trace as the test statistics. The differences in each biodiversity index in relation to the seasons and communities were assessed through linear mixed models using the pond as a random factor. The significance of the differences among seasons and among Pc, Zc and Fc was then evaluated using a Kenward-Roger adjusted F test. These analyses were performed using the functions of the "stats" (R Core Team, 2018), "lme4" (Bates, Maechler, Bolker, & Walker, 2015), "emmeans" (Lenth, 2018) and "multcompView" (Graves, Piepho, Selzer, & Dorai-Raj, 2015) packages within the R 3.4.4 programming environment (R Core Team, 2018).
Community structure of Pc, Zc and Fc in each season was evaluated using dominance-diversity curves, described by the brokenstick, niche preemption, lognormal, Zipf or Mandelbrot models of species abundance. The models imply different set of hypotheses about the niche apportionment, from lowest (preemption) to highest (lognormal and broken-stick) niche overlap, or about the relative ecological requirements of species, with species selected by environmental constraints (Zipf and Mandelbrot). Further details about the models can be found in Wilson (1991) and Magurran (2004).
The choice of the optimal model was based on Akaike's information criterion (AIC). Index calculation, curve construction and model fitting were performed using the functions of the "vegan" package (Oksanen et al., 2018) within the R 3.4.4 programming environment (R Core Team, 2018).
To estimate the functional dependence of each taxon of Pc, Zc and Fc upon the measured environmental variables, separate ridge regression models were derived for each community using the "glmnet" package (Friedman, Hastie, & Tibshirani, 2010) for the R 3.4.4 TA B L E 1 Minimum and maximum values of the variables studied during the year in the nine ponds, as well as the mean squares statistics based on the variances among ponds (MS P ) and among seasons (MS S ) Modelling of temperature-driven changes in species richness, through an approach coupling confirmatory path analysis (CPA) and non-linear modelling with artificial neural networks (ANNs), was demonstrated using environmental and D Mg data. D Mg was chosen over S ACE as the response variable in DAGs due to its lower sensitivity to parameters related to sample size and rarity of the species, upon which the latter is weighted (Magurran, 2004

| RE SULTS
All the environmental variables varied more among seasons (MS S ) than among ponds (MS P ), with the only exception of ISM (Table 1).

| D ISCUSS I ON
The  with optimal environmental conditions (Mayot et al., 2017), determined diversification of phytobenthos rather than algal blooms.
Although with similar outcomes at the ecosystem level (increase in primary productivity), the processes actually underpin very different dynamics at the community level. While algal blooms imply increases in dominance, higher slopes of the dominancediversity relationships and reductions in diversity, eventually determining mono-or oligo-specific assemblages (Gomoiu, 1992; F I G U R E 5 Directed acyclic graphs (DAGs) describing the causal relationships among T, S, W, ISM, OSM and D Mg of Pc, Zc and Fc. The AIC, Fisher's C statistic (C) and the probability (p) for each model are also reported. The model with the lowest AIC is highlighted by a black frame Green-Gavrielidis, MacKechnie, Thornber, & Gomez-Chiarri, 2018; Yoshida et al., 2018), the opposite was observed in our coastal lagoon system. Indeed, not only the number of phytobenthic taxa increased during spring, but also the dominance remained substantially the same during the year and the slope of Whittaker's curves actually decreased from winter to spring. In this context, the niche preemption model suggests a sequential build-up of the phytobenthic community, in which the number of species is limited by resource abundance rather than by interspecific competition (Magurran, 2004). This occurrence further supports the hypothesis of a resource-driven diversification of phytobenthos during spring.
It is interesting to note how limiting factors or interspecific competition (implied by the Mandelbrot community model; Wilson, 1991) shaped species abundances later, in summer, as expected considering the higher environmental constraints (water temperature, salinity, depth), eventually determining a structure typical of late successional stages in autumn, described by the broken-stick model (Magurran, 2004). A similar pattern of seasonal community development was shared by the macro-zoobenthic community, in which an analogous sigmoid model described the dominance-diversity relationship in autumn. The similarities between the dynamics of the phytobenthos and the macro-zoobenthos also encompassed the presence of two peaks of biodiversity in spring and autumn and, conversely, its decline in winter, indicating a tight dependence of the macro-zoobenthos dynamics upon those of the autotrophic group. Indeed, an increase in phytobenthos diversity would result in a greater availability of trophic and spatial niches for macro-zoobenthos and fishes (Schindler & Scheuerell, 2002;Solomon et al., 2011). However, the variability among ponds in phytobenthos diversity was widest during winter, whereas the opposite occurred in the macro-zoobenthic community, suggesting a variable versus homogeneous occurrence of the environmental constraints for the two communities across the nine coastal ponds of the "Saline di Trapani e Paceco" lagoon system. Albeit with the same seasonal dynamics in phytobenthic and macro-zoobenthic communities, both richness and diversity reached values much lower in the latter than in the former, and even lower in the fish community, indicating a decrease in biodiversity going up the trophic level. Although with some differences, this result agrees with both theoretical (Bastolla, Lässig, Manrubia, & Valleriani, 2005) and field (Cohen, Briand, & Newman, 1990;Lässig, Bastolla, Manrubia, & Valleriani, 2001) studies demonstrating biodiversity declines towards the higher trophic levels. Although seasonal variations in species richness and diversity were also observed in the fish community, with  EC 1997;Crivelli, 2006) representing the principal prey for the top predator and migrating bird Egretta garzetta in the "Saline di Trapani e Paceco" lagoon system (Piazza, Morganti, Sarà, & Campobello, 2006 Voigt et al., 2003), the same did not occur in our study system, where many phytobenthic taxa exhibited elastic net coefficients comparable or higher than those of the macro-zoobenthos and fish.
Among the candidate models relating water salinity, temperature, depth and the abundance of phytobenthos, macro-zoobenthos and fish, CPA selected a model indicating that temperature primarily affects the communities directly rather than indirectly (e.g., by changing salinity) and it is one of the main determinants of community diversity. Salinity is related to water temperature (Pawlowicz, 2013), and this variable can be also a primary shaper of marine communities (Lappalainen, Shurukhin, Alekseev, & Rinne, 2000;Williams, 1998), but does not appear to affect species richness in this system, an occurrence likely related to the abundance of euryhaline species in coastal shallow waters (Kennish & Paerl, 2010). Similarly, depth consistently affects phytobenthic taxa, mainly colonizing deeper waters, but was not selected as a driver of species richness in the three communities. Although some of the fish species observed in the lagoon system feed on zoobenthos (e.g., Pomatoschistus sp.), the model did not involve a causal relationship between species abundance of macro-zoobenthos and species abundance of fish, likely in relation to the generalist feeding behaviour of the species. The central role of temperature on species abundances in various trophic levels was already elucidated by Kirby and Beaugrand (2009) (Eviner & Chapin, 2003). Our model implies that the observed relationships between community diversity and temperature will be preserved across the +3°C ΔT above the measured temperatures. This assumption may hold for moderate warming but becomes more uncertain with more extreme increases in ΔT. Moreover, the parameterization based on one year dynamics further contributes to the uncertainty in the projections, which should be viewed as indications of possible system evolution based on its current status and as an illustration of the potential of the modelling approach in Mediterranean coastal lagoons, rather than reliable forecasts.
According to our simulations, richness of phytobenthos will show simple dynamics, with a progressive shift of the springtime peak of biodiversity to winter and the highest increase in biodiversity during summer. The diversity of the macro-zoobenthic and fish communities similarly will increase during winter and summer, an occurrence which can be related either directly to the changes in temperature or indirectly through the increase in phytobenthos diversity. Such patterns can be related to early phytobenthos development and arrival of migratory benthic or pelagic animals or their permanence for prolonged periods. Conversely, decreases in biodiversity can be easily viewed as the result of niche loss and determinants like competition and predation with newly incoming species. Indeed, biotic interactions may be the primary mediators of the impact of climate change on populations in both terrestrial (Ockendon et al., 2014) and marine ecosystems (Falkenberg, Russell, & Connell, 2012;Ghedini, Russell, & Connell, 2015;Nagelkerken, Russell, Gillanders, & Connell, 2016).
Considering the positive relationship between biodiversity and productivity in many ecosystems (Cardinale, Ives, & Inchausti, 2004;Cardinale, Palmer, & Collins, 2002;Fridley, 2001;Paine, 1966), our modelling approach may allow translating the changes to diversity across different trophic levels to changes in ecosystem productivity (Doney et al., 2012), not traditionally considered in climate change research (Tunney et al., 2014). Longer time series in the future will allow enhancing the accuracy of the model and of projections (e.g., Gauthier et al., 2013) that currently constitute a baseline and a conceptual approach to assess ecosystem structure and predict future evolutionary trajectories under climate change, allowing the development of early warning tools for conservation of entire systems.