## 2 Materials and Methods
Specifically, ELCR was calculated according to the following equation:  
\[\text{ELCR}=\text{EF}\;(\text{day}/\text{yr})\times\text{ED}\;(\text{yr}) \times\text{IR}\;(\text{kg}/\text{day})\times\text{CSF}\;(\text{mg}\;\text{kg} -1\;\text{day})\div\text{BW}\;(\text{kg})\times\text{AT} \tag{2}\]  
where EF is the exposure frequency (365 days/year); ED is the exposure duration, which was assumed to be equal to the Italian mean life expectancy (83,226 yr) (ISTAT, 2019); IR is the ingestion rate, which is equal to the PAH concentration times the mean ingestion rate of the species (FAOSTAT, 2018); CSF is the cancer slope factor for each of the analyzed PAH (OEHHA, 2009); BW is the body weight, which was assumed to be equal to the mean body weight of the Italian population (67 Kg) (ANSA, 2013); and AT is the averaging time, which is equal to EF\(\times\)ED. An ELCR above 10\({}^{-5}\), which is the acceptable lifetime risk (ALR) [52], indicates a probability greater than 1 chance over 100,000 of developing cancer [46]. THQ, which indicates the ratio between exposure and the reference dose, was calculated according to the following equation:  
\[\text{THQ}=\text{EF}\left(\text{day}/\text{yr}\right)\times\text{ED}\left( \text{yr}\right)\times\text{IR}\left(\text{kg}/\text{day}\right)\times\text{C }\left(\text{mg}/\text{kg}\right)\div\text{RFD}\left(\text{mg kg}{-1}\text{day}\right)\times\text{BW}\left(\text{kg} \right)\times\text{AT} \tag{3}\]  
where C is the PAH concentration, and RFD the oral reference dose for PAH. When THQ risk is above 1, it means that THQ is higher than the reference dose, and systemic effects may occur [46]. Graphs and statistics were produced using RStudio version3.6.1. All the analyses were carried out after grouping records according to PAHs molecular weight, i.e., analyses were conducted separately on LMW, MMW and HMW PAHs in each selected species. After checking the assumption of normality of residual distributions through the Shapiro-Wilk test (R function: shapiro_test), we used the Wilcoxon rank sum test (R function: wilcox.test) with data on Mediterranean mussel, Manila clam, red mullet and common sole to calculate the statistical significance of differences between the mean concentrations of LMW, MMW and HMW PAHs in those species. To look for seasonal trends, we used the Wilcoxon rank sum test to calculate the statistical significance of differences between the mean concentrations of each class of PAHs in cold (October-March) and warm (April-September) months. Additionally, to control for potential confounding factors, we repeated the above analysis using non-parametric ANCOVA (R function: ancova.np), with a sampling depth and sampling year as covariates. In the present work, we used the same month clustering as in [12; 53], which is based on sea water temperature: January-March (winter), April-June (spring), July-September (summer) and October-December (autumn). Finally, to see if there is a relationship between the latitude, sampling depth and sampling year and the PAH concentrations in Mediterranean mussel caught along the Italian coast of the Adriatic and the Tyrrhenian Seas, we used Kendall's rank correlation (R function: pcor; method = \"kendall\").