Species richness and vulnerability to disturbance propagation in real food webs

A central issue in ecology is understanding how complex and biodiverse food webs persist in the face of disturbance, and which structural properties affect disturbance propagation among species. However, our comprehension of assemblage mechanisms and disturbance propagation in food webs is limited by the multitude of stressors affecting ecosystems, impairing ecosystem management. By analysing directional food web components connecting species along food chains, we show that increasing species richness and constant feeding linkage density promote the establishment of predictable food web structures, in which the proportion of species co-present in one or more food chains is lower than what would be expected by chance. This reduces the intrinsic vulnerability of real food webs to disturbance propagation in comparison to random webs, and suggests that biodiversity conservation efforts should also increase the potential of ecological communities to buffer top-down and bottom-up disturbance in ecosystems. The food web patterns observed here have not been noticed before, and could also be explored in non-natural networks.

We investigated the relationships between species richness, food web complexity and the distribution of species into source and sink sub-webs. We also considered a third type, the cross sub-web, which involves disturbance propagating in both directions from an intermediate species (Fig. 1). These represent fundamental aspects of biodiversity organisation, and underlie the potential of disturbance to propagate across trophic levels 13,24,25,27,28 . Accordingly, we assumed that the higher the proportion of species co-present in at least one food chain within a food web, the higher the potential of disturbance to propagate 13 . The intrinsic vulnerability of food webs to the propagation of disturbance along food chains was thus quantified as the proportion of species (P) included in each source (P B ), cross (P C ) and sink (P T ) sub-web, originating from each basal, intermediate or top species respectively. Here, we focused on detritus-based food webs. The detritus compartment plays a fundamental role in ecosystem structure and functioning 29,30 , and climate change is expected to affect detritus inputs and organic matter decomposition rates in ecosystems 31,32 . Nevertheless, very little information that might help to understand the structure of these donor-controlled systems is available.
We compared food webs belonging to both aquatic and terrestrial habitats in the Mediterranean region (Table 1), and we quantified their vulnerability to disturbance propagating from species on various trophic levels. We then sought to verify whether the observed P B , P C and P T values were (i) dependent on species richness (S) and web connectance (Cmin) and (ii) predictable using a mechanistic food web model 6 or similar to what would be expected by chance. This was achieved by comparing field data with random food web models, as well as webs generated using the niche model 6 , a simple yet predictive food web model based on observed S and Cmin values and a few foraging rules constraining the probability of feeding interaction between species. The comparison of observed and random food web structures can provide important insight into the organisation of ecological communities. Indeed, the association of food web stability with non-random structures enhances our understanding of the ecological factors constraining the observed patterns and hence the mechanisms underlying the emergence and persistence of complex food webs 7 . and trophic links (lines). In network 1, all species are part of one food chain, and a disturbance originating from any given species would directly propagate along the food chain to all the remaining ones. In network 2, a bottom-up disturbance spreading from species 1 would affect species 4-5-6 only. In network 3, the modification of a single link with respect to network 2 would mean that only species 4 would be directly affected by a bottom-up disturbance propagating from species 1. Panels (b-d) Food sub-webs exemplifying various propagation pathways along food chains for disturbance starting from a single species (circled): (b) bottom-up propagation from a basal resource (red), (c) cross propagation from a primary consumer (orange), and (d) topdown propagation from a predator (yellow).

Results
Species richness and food sub-webs. Food web connectance (Cmin) varied from 0.073 in Stream to 0.130 in Lake (Table 1), decreasing with species richness (S) (Fig. 2). On average, bottom-up propagation pathways (P B ) included a higher proportion of species than cross (P C ) and top-down (P T ) propagation pathways ( Fig. 3; one-way ANOVA for repeated measures: F = 36.7, p < 0.0001; associated Tukey's pairwise comparisons: Q always > 4.4, p always < 0.01). P B , P C , P T and V values differed between habitats (Fig. 3), all of them increasing with Cmin, while decreasing with S and intermediate species as a percentage of the total (%I) ( Table 2 and Fig. 4). In contrast, P B , P C and P T were not related to the proportion of basal or top species (p always > 0.05). No differences in the correlation slopes between P B , P C and P T and either S or %I were observed (ANCOVA and the  Table 1. Food webs representative of various habitats in the Mediterranean region. Locality indicates the region, country and GPS coordinates for each habitat. S: number of species in the food web; Cmin: food web connectance; L/S: species' feeding linkage density.

Figure 2.
Panels (a,b) relationship between species richness (S), number of feeding links (L) and food web connectance (Cmin) observed in our dataset (red symbols) and estimated from the link-species scaling law (black symbols) and the constant connectance hypothesis (empty symbols). Panels (c,d) observed mean P B , P C and P T values (red symbols) in comparison with values estimated from the link-species scaling law (black diamonds) and the constant connectance hypothesis (empty diamonds), as well as values measured in model food webs generated using the niche model (black triangles) and a random model (empty triangles). P represents the proportion of species in each source (P B ), cross (P C ) or sink (P T ) food sub-web. Lower panels (e-h) relationship between values observed in real food webs and those measured in niche model-generated webs. V represents the mean Pi value of all nodes in the food web. Grey lines represent y = x ± 0.15. Given that the maximum possible value of P is 1, the grey lines represent a ± 15% deviation from y = x. (2019) 9:19331 | https://doi.org/10.1038/s41598-019-55960-8 www.nature.com/scientificreports www.nature.com/scientificreports/ homogeneity of slopes test, F < 1 and p > 0.05 in both cases). In contrast, the correlation slope between P B and Cmin was lower than what was observed for P C and P T (ANCOVA and homogeneity of slopes test, F = 4.5 and p = 0.01).

Figure 3.
Comparison of the proportion of species (P) in the three types of disturbance propagation pathway. P B , P C and P T represent the mean proportion of species belonging to source, cross and sink sub-webs respectively. Food webs are ranked in order of V value, which represents the mean P value of all species within each food web.  Table 2. Multiple regression models testing the effect of number of species (S), intermediate species as a percentage of the total (%I) and food web connectance (Cmin) on food web vulnerability to bottom-up (P B ), cross (P C ) and top-down (P T ) propagation of disturbance across trophic levels, as well as on overall vulnerability to disturbance propagation (V). Bold values indicate a significant effect, with a p value < 0.05.

Discussion
The distribution of species across source, cross and sink sub-webs differed between habitats in a predictable way in accordance with food web size (S) and connectance (Cmin). The higher the number of species, the lower the proportion of them co-present in at least one food chain. This indicates that the vulnerability of ecosystems to disturbance propagation across trophic levels falls as their species richness rises, which implies that more biodiverse communities are potentially more able to buffer disturbance. Our observations are also consistent with recent results from studies based on model food webs of similar species richness 13 , where higher connectance intensified the effect (i.e. caused a higher rate of secondary species loss) of single perturbation events (i.e. successful species invasions). As well as species richness, a high number of intermediate species relative to the total reduced the potential of disturbance to propagate along food chains. In addition, our data satisfied expectations from the link-species scaling law, which predicts constant linkage density among species. According to foraging optimisation theories, (i) a low and roughly constant linkage density between species arises from the selection by consumers of the lowest number of food items that maximises their net energy intake 9,35 , and (ii) prey species richness promotes predator trophic specialisation, as shown in the food webs analysed here 32,36 and elsewhere 9,37 . This suggests that energetic constrains on foraging strategies, which operate at the individual level 38 , may give rise to food webs which are intrinsically less vulnerable to the propagation of disturbance than what would be expected if consumers tended to generalise as the number of their potential food sources increased 39 . Notably, the random food web models failed to reproduce real patterns, while niche model-generated webs closely predicted the observed values. Our results thus imply that the distribution of species in source, cross and sink sub-webs significantly deviates from random, and that the patterns observed in real food webs benefit ecological communities by reducing vulnerability to disturbance propagation along food chains compared to what would be expected if distribution was determined by chance. Regardless of habitat type, source sub-webs included a higher proportion of total species than cross and sink sub-webs. This implies that disturbance propagating from a basal resource would have a significantly higher potential to spread throughout the food web than disturbance propagating from a top or an intermediate species.
We speculate that the observed patterns may promote stability in detritus-based systems, where the availability of basal resources (i.e. organic detritus and colonising microfungi) is expected to be much more stable over time than invertebrate populations, often characterised by marked spatial-temporal variations and sensitivity to stress 40,41 (see Fig. S2 for supporting results). On the other hand, global change is expected to affect the dynamics of the detritus compartment in ecosystems 31,32,42,43 . Thus, our data suggest that there is a strong likelihood of future environmental changes in Mediterranean ecosystems giving rise to bottom-up effects mediated by the structure of detritus-based food webs.
The coefficient of variation of Pi values increased with the number of species within each food web, and it was inversely related to food web vulnerability to disturbance. Together with the P B , P C , P T and V values, this makes the variation of species' Pi values within a food web a useful a-priori indication of its vulnerability. We acknowledge that our observations are based on static food web structures, and that the reorganisation of trophic links and/or changes in linkage strength between species in accordance with their trophic niche plasticity could www.nature.com/scientificreports www.nature.com/scientificreports/ modify the effects of disturbance 25 . Nevertheless, our analysis of sub-webs enabled us to evaluate the intrinsic vulnerability of real food webs to a range of disturbance propagation pathways (i.e. bottom-up, cross, and top-down). This approach is expected to yield information useful to management strategies based on the risk of ecosystem-specific disturbance propagation, regardless of food web compilation methods. Analysis of the scaling relationships between the size of networks and that of the sub-webs of which they are composed could also be extended to non-natural systems (e.g. financial 44 , transport 45 or disease vector dispersal 46 networks). This would make it possible to verify whether (i) systemic risk increases with network size or (ii) the rules governing such systems mean that an increase in size is accompanied by increased potential to buffer disturbance, as observed in the food webs analysed in this study.
In conclusion, we show that increasing species richness and foraging constraints on food web complexity limit the proportion of species that share food chains in ecosystems. This pattern has not been noticed before. It reduces the intrinsic vulnerability of biodiverse ecological communities to disturbance propagation, and promotes the emergence of species-rich food webs where species organise into effective disturbance-buffering structures in accordance with predictable rules. It is also consistent with the predictions of theoretical research, which highlight the positive effect of adaptive trophic behaviour by consumers on the stability of model food web networks 9 . Conversely, the observed results imply that species-poor ecosystems could be highly vulnerable to the propagation of disturbance along food chains, including bottom-up disturbance associated with climate change 31,42 , with increased risks for the stability of food webs and the ecosystem services they support [47][48][49] . While this study focused on detritus-based food webs, the observed results may be extended to more complete food webs including primary producer-based food chains. Indeed, detrital and herbivore energy pathways are closely interconnected 8 , and our observations are consistent with webs generated using the niche model, which was developed and tested against whole food web structures in both aquatic and terrestrial ecosystems 6 . Notably, our food web analysis provides mechanism-based evidence indicating that efforts devoted to biodiversity conservation also increase the potential of natural communities to buffer disturbance in ecosystems by maintaining biodiverse, relatively less complex, and 'safer' food web structures for their constituent species.

Materials and Methods
We compared habitat-scale detritus-based food webs (Table 1), reconstructed using population abundance data, carbon and nitrogen stable isotope values and Bayesian isotopic mixing models (for details on stable isotope analyses and food web reconstruction methods, see Rossi et al. 32 and Calizza et al. 12,36 ). Specifically, the habitat-scale food webs were based on taxa and trophic links obtained by sampling multiple patches within each habitat. This made it possible to taking into account intra-habitat variability of occurrence data and the trophic preferences of the various taxa 12 . The food webs were representative of six different habitats: a river flowing from a semi-natural to an anthropic landscape (the River Sacco, referred to as Stream), a lowland high-order river embedded in an urban landscape crossing the city of Rome (the River Tiber, referred to as River), a meso-oligotrophic lake (Lake Bracciano, referred to as Lake); a transitional water system (the Santa Gilla lagoon, referred to as Lagoon), arable land (referred to as Corn Field) and a woodland (referred to as Beech Forest) ( Table 1). The lower stretches of the two river courses were more heavily anthropised and barriers were present in both cases (e.g. the city of Rome along the River Tiber). Therefore, River and Stream were each divided into upstream and downstream habitats. In all habitats, basal resources were represented by habitat-specific leaf litter colonised by various species of microfungi 26,50 . Detailed descriptions of sampling activities and habitats can be found in Rossi et al. 32 (Lake, Lagoon, Beech Forest, Corn Field), in Calizza et al. 12,36 (River), and in Bentivoglio et al. 51 (Stream).
Web topology was analysed by determining (1) species richness, S; (2) the number of feeding links, L; (3) linkage density, L/S; (4) connectance (Cmin) as L/S 2 ; and (5) the proportions of total species in the food web accounted for by basal (%B), intermediate (%I) and top (%T) species. %I was considered to be a measure of prey availability for predators 32,36 . The intrinsic vulnerability of food webs to bottom-up, cross and top-down propagation of disturbance (i.e. starting from basal, intermediate and top species) was quantified as the proportion of the total number of species (P) occurring in each source (P B ), cross (P C ) and sink (P T ) sub-web respectively 24 (Fig. 1), as follows: where nx is the number of basal, intermediate and top species used to calculate P B , P C and P T respectively, p i = Si/ (S-Ni), Si = the number of species included in the sub-web originating from the i th species, S = the total number of species in the food web, and Ni = number of species on the same trophic level as the i th species. The overall vulnerability of the food web to disturbance propagation along food chains (V) was measured as the mean P value of all species in the food web.
For each parameter, pairwise ratios between habitats were calculated and a matrix containing pairwise comparisons was created. Pairwise values (n = 28) were plotted and used for correlation models 52,53 . A permutation test (9999 permutations) was also run on correlation coefficients, and the observed correlations were considered significant only if robust to permutation (i.e. with a permutation-based p value < 0.05). Differences in correlation slopes were tested by analysis of covariance (ANCOVA) and the associated homogeneity of slopes test.
Based on the observed S and Cmin values, for each habitat, 20 model food webs were generated. Specifically, 10 model food webs were generated based on the niche model 6 , a simple yet predictive one-dimensional model able to reproduce the topological properties of real food webs starting from a small number of rules governing the probability of feeding interaction between species. The other 10 model food webs were generated by assigning links between species at random. The P B , P C , P T and V values were calculated and compared with the values measured in the real food webs. In addition, the L and Cmin values expected in accordance with the link-species