*2.4. Data Analysis*

Di fferences between lakes in terms of animal community composition (considering both fish and benthic invertebrates) were tested using contingency tables based on chi-square ( χ2) tests, Monte Carlo permutation tests and the associated Cramer's V index (a measure of the strength of association among communities; Past 3.0 software package). Specimens collected in each sampling site (and replicates) were grouped by type or taxon (respectively for basal resources and animals) for each lake.

The Shannon diversity index (Hs) of invertebrate fauna for each lake was calculated at family level considering a total abundance of the taxa collected in each lake. Given that assessing Shannon diversity is only possible at the level of equal identification of all taxa, the few individuals belonging to the Gastropoda, Oligochaeta, Nematoda, and Nemertea classes (together accounting for less than 1.5% of total fauna) were excluded from the Shannon diversity index computation. Hutcheson's diversity *t*-test and the associated bootstrap procedure (9999 replicates), both available in the Past 3.0 software package, were applied to Hs values to test for significant di fferences [51]. Hutcheson's diversity *t*-test is a modified version of the classic *t*-test and is based on comparison of Hs variances. The *t* statistics of Hutcheson are defined as: 

$$t = \frac{\left|Hs\_i - Hs\_j\right|}{\sqrt{var(Hs\_i) + var\left(Hs\_j\right)}}\tag{1}$$

which follows Student's *t* distribution. In the equation, *i* and *j* referred to the invertebrate communities of the lakes in paired comparisons, *Hs* represents the Shannon diversity index and *var(Hs)* its variance.

The isotopic values of collected organisms were used to reconstruct the diets of the eel *Anguilla anguilla* and the annular seabream *Diplodus annularis* in each lake. The diets were estimated on the basis of the Isotopic Trophic Unit (ITU) method [31]. The isotopic signatures of single basal resources, invertebrates, and fish were represented in the bi-dimensional isotopic space (Figure S1). This was subdivided into squares (ITUs) corresponding to 1 × 1% δ15N and δ13C values, starting from the lowest δ13C value in the dataset and a δ15N value of zero. The ITUs were thus identified and labelled (Figure S1).

The diets of each ITU containing individuals of the two fish species were calculated by means of Bayesian Mixing Models (R software ver. 3.5.3, SIMMr package) [52] considering a Trophic Enrichment factor (TEF) of 3.4 ± 1.0% for δ15N and 1.0 ± 0.5% for δ13C [18,37,49,53–56] and uninformative priors. These TEF values (expressed as mean ± standard deviation) are considered a robust and widely applicable assumption in the presence of multiple trophic pathways and di fferent types of food sources [37,56]. For all SIMMr models, we ran three Markov Chain Monte Carlo chains of 300,000 iterations each with a burn-in of 200,000 and a thinning rate of 100 iterations. We assumed that all incoming food items had the same probability of being included in the consumer's diet. The model considers both variance in the isotopic signatures of the resources and uncertainty regarding the trophic enrichment of the consumer (TEF). The model results were expressed in the form of a probability distribution of plausible contribution values. The central tendency values of the distribution (mode, mean, median) allowed us to identify the most important food sources, while the upper and lower limits of the credibility ranges (CI: 50%, 75%, 95%) revealed the range of feasible contributions. The pool of food sources was selected based on the mixing model outputs in accordance with Rossi et al. [31]. Since the *A. anguilla* and *D. annularis* diet was obtained by starting from the foraging choice of each individual, the overall contribution of some food sources, important at individual level (>5%), could be relatively small (<5%) if considered at the population level (see also [31]). In order to obtain detailed information on the diet of the eel and the seabream these contributions were also considered.

Individuals other than *A. anguilla* and *D. annularis*, including basal resources and invertebrates, were excluded from ITU-consumers (but not from potential ITU food sources) before performing the Bayesian mixing models. This was done in order to correctly estimate the diet of *A. anguilla* and *D. annularis*. The set of potential ITU food sources was considered on the entire δ13C axis and within a given range on the δ15N axis, i.e., within ±3.4% (the TEF) of the value of the consumer [31]. The Bray–Curtis similarity index (BC), based on the contribution of each resource to the diet of the two fish, was also calculated in order to quantify the diet similarity among lakes [36,56]. BC is expressed as proportional similarity ranging from 0, when no common food sources are found for the compared groups, to 1, when the compared groups have the same food sources in the same proportions [36,56].

The symmetric overlap in resource use [57–59] was measured in accordance with the Pianka equation [59]:

$$O\_{jk} = \frac{\sum\_{i=1}^{n} p\_{ij} p\_{ik}}{\sqrt{\sum\_{i=1}^{n} \left(p\_{ij}\right)^2 \sum\_{i=1}^{n} \left(p\_{ik}\right)^2}} \tag{2}$$

where the Pianka index ( *Ojk*) represents a symmetric measure of overlap between species *j* and *k*, and *pij* and *pik* are the proportional contributions of any given resource *i* used by species *j* and species *k*. The Pianka index ranges from 0 (overlap absent) to 1 (complete overlap).

Chesson's selectivity index [60] was calculated for each food item to determine possible preferences for particular food sources among those o ffered:

$$\alpha\_{i} = \frac{r\_{i}/n\_{i}}{\sum\_{i=1}^{m} (r\_{i}/n\_{i})} \tag{3}$$

where αi is Chesson's selectivity index, *m* is the number of food source types, *ri* is the proportion of food type *i* in the diet and *ni* is the proportion of food type *i* in the environment. The value of α*i* ranges from 0 to 1, with 0 indicating complete avoidance, values above 1/m indicating preference and 1 indicating absolute preference [61]. Since consumer isotopic ratios provide an integrated measure of prey assimilated over time, we hypothesized that the composition of the taxon in each lake did not vary considerably over the course of a season. Therefore, the Chesson index based on the relationship between assimilated prey and its abundance in the environment could measure the selectivity of food products with a good approximation.

χ2 tests were performed to test for differences between lakes in terms of the relative abundance of fauna and differences between food sources in terms of their proportional contribution to the diet of the fish population in each lake. Although it is not possible to establish a theoretical expected value, a χ2 test was performed to test for possible differences between IP and LP and between HP and LP, considering the least polluted lake as the reference value.

Differences between lakes in both the δ15N and δ13C isotopic signatures of basal resources and fauna were tested by one-way ANOVA for comparisons between normal distributions (Shapiro–Wilk normality test, *p* > 0.05) while the Mann–Whitney with Bonferroni correction in cases of multiple comparisons and Kruskal–Wallis tests were used if non-normality was observed (Shapiro–Wilk normality test, *p*-value < 0.05). Levene's test for variances was used to test for differences within and between lakes in the δ13C variance of primary producers. Kruskal–Wallis tests and associated Mann–Whitney pairwise comparisons were also used to compare the proportional contribution of food items to the diet of *A. anguilla* and *D. annularis*.

The niche metrics for both species in each lake were also calculated [62–64]. These metrics, originally proposed by Layman et al. [62] for application at community level, can be used at population level to obtain information about trophic diversity within a single population [35,63,64]. These included the ranges (highest to lowest) of δ13C (Carbon Range, CR) and δ15N (Nitrogen Range, NR) values. CR provides information about the variety of food sources exploited by the population (i.e., its trophic generalism), while NR indicates the number of trophic levels (i.e., degree of omnivory) of the population. The isotopic niche widths of both *A. Anguilla* and *D. annularis* were calculated as SEAc (Standard Ellipse Area corrected by degree of freedom) using R software ver. 3.5.3, SIBER analysis package [64,65]. The SEAc encompasses the core (about 40%) of the population's isotopic observations. This is a solid metric for comparing the isotopic niche of populations regardless of sample size and any isotopic outliers in the data [62,64]. Linkage density (L/S) was measured as the average number of feeding links (L) per ITU (S). Finally, based on the proportional contribution of each food source, the trophic niche width (TNW) of each population was measured as the diversity of resources consumed (Hs) by each population and compared among lakes. If not specified otherwise, the results are reported as mean ± standard error (s.e.).
