Safety Evaluation and Probabilistic Health Risk Assessment of Cow Milk Produced in Northern Italy According to Dioxins and PCBs Contamination Levels

Contamination levels of dioxins and polychlorinated biphenyls (PCBs) were monitored over 2018–2021 in 214 bovine milk samples from farms located in two regions in northern Italy (Lombardy and Emilia-Romagna). The average concentrations of the sum of dioxins and dioxin-like PCBs (0.78 ± 0.55 pg TEQ/g fat) and six non-dioxin-like PCBs (6.55 ± 2.24 ng/g fat) were largely below the maximum, and action limits established at European level, confirming a decreasing trend observed both locally and across Europe in recent years. The impact of contamination levels on chronic dietary exposure of the Italian population to dioxins and PCBs was found to be highly variable based on the type of cow milk (skimmed, semi-skimmed, or whole-fat milk) and the population age group considered. Indeed, a first-tier screening of the potential exposure via determinist methods allowed for the identification of the youngest population as the group with the worst risk profile. The refinement of exposure assessment via Monte Carlo probabilistic methods suggested that, at the less pessimistic middle-bound simulation scenario, infants, toddlers, and children consuming whole cow milk may be exposed to dioxins and PCBs levels above the toxicological reference values with a probability of 76, 56, and 22%, respectively.


Introduction
Polychlorinated dibenzo-p-dioxins and furans (PCDD/Fs, also known as dioxins) and polychlorinated biphenyls (PCBs) are persistent organic pollutants widely distributed in the environment. Even though they include a wide range of chemical compounds, they all have low degradation potentials and high lipophilic character. Furthermore, PCDD/Fs and PCBs have similar chronic toxic effects on animals and humans, potentially causing developmental and reproductive disorders, endocrine and immune system impairment, teratogenesis, and tumorigenesis [1][2][3][4]. lowered the MLs for several food commodities in 2022. For milk and dairy, MLs of 2 and 4 pg TEQ/g fat of PCDD/Fs and PCDD/Fs + dl-PCB, respectively, and of 40 ng/g fat for non-dioxin-like (ndl)-PCBs are currently in force [32].
In light of new TWI and MLs limits, as well as changes in the dietary habits of the Italian population, it is therefore important to pursue risk management strategies through the ongoing monitoring of dioxin and PCB contamination levels in food and feed to measure the public health impact of dietary exposure to these contaminants and to identify those segments of the population that are at increased risk of exposure and potentially long-term toxic effects. In this context, probabilistic risk assessment might provide more confident and realistic representations of the actual risk levels to which consumers are exposed, accounting for the variability and uncertainty of exposure estimates.
In this work, the results of the 2018-2021 monitoring plan on dioxin and PCB residues carried out in the northern Italian regions of Lombardy and Emilia Romagna are presented, with the aim of evaluating the occurrence, spatial distributions, and temporal variations of these contaminants in bovine milk produced in this territory. A second-level objective was to perform an exposure assessment of dioxins and PCBs intake through the consumption of milk with different fat contents by different population age groups in Italy and verify to what extent the consumption of milk contributes to the toxicological reference values using both deterministic and probabilistic (Monte Carlo method) approaches, providing a realistic scenario of the possible risks.

Milk Samples
Milk samples were collected from 2018 to 2021 under the framework of the Italian residue monitoring plan (RMP). In those years, the Italian RMP was structured according to Council Directive 96/23/EC [33] and Commission Decision 97/747/EC [34], as well as Regulation (EU) 625/2017/EU [35]. A total of 214 bovine milk samples were obtained from 187 dairy farms in the Lombardy and Emilia-Romagna regions of Northern Italy (Figure 1). Milk samples were collected from bulk farm tanks by the official competent authorities according to Regulation (EU) 644/2017/EU [36]. Specifically, bulk milk was thoroughly mixed in order to ensure and homogeneous distribution of the contaminants, and an aggregate sample of 1 L was obtained from each farm. Once collected, the samples were stored at +4 • C and delivered within 24
Pre-packed multi-layer silica, alumina, and carbon columns were produced by FMS (Fluid Management System, Billerica, MA, USA). The 13 C-labeled recovery, clean-up, and standard injection solutions were provided by CIL (Cambridge Isotope Laboratories, Andover, MA, USA). For PCDD/Fs, EDF-9999 Method 1613 calibration solutions (CS1-CS5) were used. For PCB calibration, an in-house curve was prepared using PCB MIX-75
All solvents and reagents used for the analyses were tested to ensure the absence of contaminants at the levels of interest (i.e., below one-fifth of the limit of quantification (LOQ) for PCDD/Fs and below one-tenth of the LOQ for dl-PCB and ndl-PCB).

Analytical Method
The methods used for the determination of PCDD/Fs and PCBs in milk samples were previously described by Lorenzi et al. [20,37]. Briefly, milk samples were homogenized and freeze-dried (Freeze Dryer Martin Christ Gefriertrocknungsanlagen, Osterode am Harz, Germany). Following further homogenization, a portion of each sample was tested with the Soxhlet method to determine the lipid content, while 8-10 g of the samples were mixed with diatomaceous earth and spiked with the 13 C-labeled internal standard EDF-8999 (Cambridge Isotope Laboratories, Andover, MA, USA) and the 13 C-labeled internal standard EC-4995 (Cambridge Isotope Laboratories, Andover, MA, USA). An accelerated solvent extractor (ASE ® 300, Sunnyvale, Dionex, CA, USA) was used for fat extraction with toluene (2 cycles of 20 min at +135 °C and 1500 psi). The solvent was filtered using anhydrous sodium sulfate and evaporated at +50 °C with a rotatory evaporator. The lipid extracts were dried overnight at 70 °C ± 5 °C. Fat content was determined gravimetrically and compared to the value obtained from the Soxhlet method to confirm ASE extraction efficiency.
Lipid extracts were solubilized in 5 mL of hexane/dichloromethane solution (1:1, v/v), spiked with the standard clean-up solution EC-4978 (Cambridge Isotope Laboratories, All solvents and reagents used for the analyses were tested to ensure the absence of contaminants at the levels of interest (i.e., below one-fifth of the limit of quantification (LOQ) for PCDD/Fs and below one-tenth of the LOQ for dl-PCB and ndl-PCB).

Analytical Method
The methods used for the determination of PCDD/Fs and PCBs in milk samples were previously described by Lorenzi et al. [20,37]. Briefly, milk samples were homogenized and freeze-dried (Freeze Dryer Martin Christ Gefriertrocknungsanlagen, Osterode am Harz, Germany). Following further homogenization, a portion of each sample was tested with the Soxhlet method to determine the lipid content, while 8-10 g of the samples were mixed with diatomaceous earth and spiked with the 13 C-labeled internal standard EDF-8999 (Cambridge Isotope Laboratories, Andover, MA, USA) and the 13 C-labeled internal standard EC-4995 (Cambridge Isotope Laboratories, Andover, MA, USA). An accelerated solvent extractor (ASE ® 300, Sunnyvale, Dionex, CA, USA) was used for fat extraction with toluene (2 cycles of 20 min at +135 • C and 1500 psi). The solvent was filtered using anhydrous sodium sulfate and evaporated at +50 • C with a rotatory evaporator. The lipid extracts were dried overnight at 70 • C ± 5 • C. Fat content was determined gravimetrically and compared to the value obtained from the Soxhlet method to confirm ASE extraction efficiency.
Lipid extracts were solubilized in 5 mL of hexane/dichloromethane solution (1:1, v/v), spiked with the standard clean-up solution EC-4978 (Cambridge Isotope Laboratories, Andover, MA, USA), containing three 13 C-labeled PCB congeners, and diluted with 20 mL of hexane. The diluted extracts were purified using pre-packed silica columns, acidified with sulfuric acid, and eluted with n-hexane. The purification fractions were concentrated to 0.5 mL in the TurboVap evaporator (Zymark, Hopkinton, MA, USA) and loaded into the Power-Prep ™ system (Fluid Management System, Watertown, MA, USA) equipped with silica, alumina, and carbon columns. Toluene (50 mL) was used to elute PCDD/Fs from the carbon column, while n-hexane and a mixture of hexane/dichloromethane solution (9:1, v/v) were used for PCB elution from the alumina column. The final extracts were dried using the TurboVap evaporator and a vacuum concentrator (Genevac, Ipswich, UK).
The PCDD/F fraction was dissolved in 10 µL of 1:100 ED-2521 injection solution (Cambridge Isotope Laboratories, Andover, MA, USA), and the PCB fraction was dissolved in 20 µL of 1:50 EC-4979 injection solution (Cambridge Isotope Laboratories, Andover, MA, USA). HRGC-HRMS analysis was carried out using two TRACE GC ULTRA gas chromatography coupled with a DFS (double focusing system) high-resolution magnetic scan system (Thermo Fisher Scientific, Waltham, MA, USA) or with AutoSpec high-resolution mass spectrometer (Micromass/Waters, Manchester, UK) [24]. A DB5 MS capillary column (60 m × 0.25 mm, 0.25 µm, Thermo Fisher Scientific, Waltham, MA, USA) was used for PCDD/F separation, and a TR-PCB 8MS capillary column (50 m × 0.25 mm, 0.25 µm, Thermo Fisher Scientific, Waltham, MA, USA) for PCB separation. The mass spectrometer operated in selected ion monitoring mode with a mass resolution of over 10,000.
Results were expressed in pg/g fat for PCDD/Fs and dl-PCBs and ng/g fat for ndl-PCBs. Toxic equivalent values (TEQ) were calculated using the 2005 World Health Organization Toxic Equivalency Factors (WHO-TEFs) [5].

Quality Control
The laboratory participated in different international proficiency tests and interlaboratory studies, obtaining satisfactory results.
To ensure the absence of cross-contamination and interferences, blank samples were processed with each batch of five samples. For every ten samples, a fortified sample was processed to ensure analyte recovery. The stability of response factors for all congeners was checked daily using calibration verification standard solutions for PCDD/Fs and PCBs and assessing the acceptability of deviations from the calibration curves [38,39]. The recovery values of the internal standard congeners were determined for each sample.
Duplicate analyses were performed on each sample, and the compliance of the absolute difference between the two measurements with the repeatability limits, obtained during the validation study, was used to estimate the overall repeatability of the method.
LOQ and recovery values are reported in Supplementary Table S1.

Summary and Statistics
Contamination levels of dioxins and dl-PCBs were reported both as actual mass fraction values (pg/g fat) and toxic equivalency values (pg TEQ/g fat), while those of ndl-PCBs as actual mass fraction values (ng/g fat).
Non-detected values, i.e., values below the LOQ, were substituted with the LOQ (upper bound concentrations, UB), 1 /2 LOQ (middle bound concentrations, MB), or 0 (lower bound concentrations, LB). Univariate descriptive statistics were applied to the UB-data matrix to identify potential differences among the tested contaminants. The Student's T test was employed to evaluate differences between the two Italian regions from which samples were collected (Lombardy and Emilia-Romagna). The Analysis of Variance (ANOVA) followed by Tukey's post hoc test was used to identify the provinces of each Italian region characterized by the highest concentrations of contaminants and to explore significant differences over the 2018-2021 sampling years. A Spearman's correlation analysis was performed to investigate potential links between the fat content and dioxin concentrations in milk samples. PCDD/F + dl-PCB concentrations (pg/g fat) of each sample were multiplied by the relative fat content (g fat/100 g) to calculate the overall PCDD/Fs + dl-PCBs in 100 g of milk, and then these results were correlated with the respective lipid content. Correlation coefficients (r) higher or lower than 0.6 were considered indicative of the existence of a positive or negative correlation, respectively. The statistical significance was set at p ≤ 0.05 for all tests.

Screening via Deterministic Approach
The estimated weekly intakes (EWIs) of PCDD/Fs + dl-PCBs through milk consumption were calculated using the measured PCDD/Fs + dl-PCBs mean (UB) concentrations, milkfat concentrations, and mean or 95th percentile chronic dietary consumption data of cow milk by the Italian population.
Consumption data of cow milk (consumers only) were retrieved from a more recent Italian survey on food consumption SCAI IV included within the EFSA Comprehensive Food Consumption Database [40][41][42]. Infants (up to 12 months old), toddlers (13-36 months old), children (37 months-9 years old), adolescents (10-17 years old), adults (18-65 years old), and the elderly (>65 years old) were the six Italian population age groups considered.
The cow milk collected in this study was supposed to be destined for the production and commercialization of milk with different fat percentages. Therefore, dietary consumption data of skimmed, semi-skimmed, and whole cow milk were retrieved from Level 6 of the FoodEx2 hierarchical classification system included in the EFSA Database and used to simulate three different chronic exposure scenarios per population group. For the calculation of EWIs, skimmed, semi-skimmed, and whole cow milk were supposed to contain 0.3, 1.8, and 3.5 g lipids/100 g, respectively, by using Equation (1): where C is the mean concentration of PCDD/Fs + dl-PCBs (pg TEQ/g fat); M is the mean or 95th percentile chronic weekly consumption of skimmed, semi-skimmed, or whole cow milk (g/kg bw/week); F is the fat content of skimmed, semi-skimmed, or whole cow milk samples (g fat/100 g). Finally, the risk was characterized by comparing the calculated EWI with the TWI of PCDD/Fs + dl-PCBs (equal to 2 pg/kg bw/week) by using Equation (2): (2)

Refinement via Probabilistic Methods (Monte Carlo Simulations)
The Monte Carlo (MC) probabilistic method was applied to simulate the distribution probabilities of EWI values and the related %TWI only for those deterministic scenarios indicating potential causes of concern, as suggested by EFSA [43]. Rather than using single point estimates, the distribution frequency of both the PCDD/Fs + dl-PCB and chronic consumption data of cow milk (i.e., assumptions) were considered. In this way, information concerning uncertainty, variability, and probability in health risk evaluation was obtained.
First, the best-fit probability distributions for the PCDD/Fs + dl-PCB contamination levels were tested using the modified Kolmogorov-Smirnov goodness-of-fit test (p > 0.05). In this case, the LB, MB, and UB concentrations were individually used to account for uncertainty embedded within the data. Based on the results, the LB, MB, and UB concentrations were fitted to the lognormal distribution. No distribution frequencies were set for milkfat, and fixed values at 1.8 and 3.5 g/100 g for semi-skimmed and whole cow milk were used. The overall distribution frequency of dietary consumption data for milk (expressed as g/kg bw/week) was supposed to be lognormal based on the literature data [44,45].
After that, MCSs were run randomly and repeatedly using 50,000 iterations for each population subgroup using Equations (1) and (2) so as to forecast a set of theoretically stable distributions of PCDD/F + dl-PCB EWIs and relative contributions to TWI (%), as well as their probabilities at different risk levels (10-90th percentiles).

Spatial and Temporal Distribution of Dioxins and PCBs in Milk
The average concentrations of dioxins and PCBs measured in the 214 cow milk samples are reported in Table 1, where they were statistically summarized by the Italian region (i.e., Lombardy vs. Emilia-Romagna) and the year of sampling (i.e., 2018 vs. 2019 vs. 2020 vs. 2021). Individual concentrations of each of the 35 measured PCDD/F, dl-PCB, and ndl-PCB congeners are provided in Supplementary Tables S2 and S3. Regardless of both the region and the year, the average concentrations of PCDD/Fs, dl-PCBs, PCDD/Fs + dl-PCBs, and ndl-PCBs were found to be 0.21 pg TEQ/g fat, 0.54 pg TEQ/g fat, 0.75 pg TEQ/g fat, and 0.64 ng/g fat, respectively (Table 1).   No statistically significant differences (p > 0.05) concerning concentrations in all the tested contaminants were found according to the four years of sampling. On the contrary, milk samples from Emilia-Romagna were found to be significantly less contaminated (p ≤ 0.05) by both dioxins and PCBs than those collected from the Lombardy region. Indeed, as can be observed in Table 1, milk samples from Emilia-Romagna presented a 20% lower concentration of PCDD/Fs (0.19 vs. 0.23 pg TEQ/g fat), a 40% lower concentration of dl-PCBs (0.41 vs. 0.66 pg TEQ/g fat), 30% lower concentrations of PCDD/Fs + dl-PCBs (0.60 vs. 0.88 pg TEQ/g fat), and 7% lower concentrations of ndl-PCBs (6.25 vs. 6.71 ng/g fat).
A deeper analysis of data was conducted to investigate whether the increased contamination of milk with PCDD/Fs + dl-PCBs was proportional to the lipid content of milk. A Spearman correlation analysis between the sample contamination levels and the lipid content was performed. Interestingly, the analysis showed a positive correlation between contamination levels and fat content only in milk samples from Emilia-Romagna (r = 0.63, p ≤ 0.05), while no significant correlations were found in samples from Lombardy (r = 0.11, p > 0.05). These findings suggest that Emilia-Romagna may have a relatively consistent and stable level of background contamination to the extent that the fat content of milk from this region could be somewhat indicative of the degree of dioxin contamination. In contrast, milk samples from Lombardy may be affected by background and point contamination, which could explain why some samples with low-fat content showed high contamination levels while some samples with high-fat content showed low contamination levels. This situation may be attributed to the higher level of industrialization in Lombardy, the existence of contaminated sites that produced PCBs in the past, and the proximity of industrial facilities to dairy farms [46]. All milk samples were compliant with the European MLs of PCDD/Fs of 2 pg TEQ/g fat and ALs of 1.75 pg TEQ/g fat, as well as with the European MLs of ndl-PCBs of 40 ng/g fat. Nonetheless, one sample from Lombardy exceeded the MLs of 4 pg TEQ/g fat established for PCDD/Fs + dl-PCBs since contaminated by 6.28 pg TEQ/g fat, while the other four samples from Lombardy were above the ALs of 2 pg TEQ/g fat established for dl-PCBs ( Figure 2). from this region could be somewhat indicative of the degree of dioxin contaminati contrast, milk samples from Lombardy may be affected by background and point co ination, which could explain why some samples with low-fat content showed high tamination levels while some samples with high-fat content showed low contamin levels. This situation may be attributed to the higher level of industrialization in bardy, the existence of contaminated sites that produced PCBs in the past, and the imity of industrial facilities to dairy farms [46].
All milk samples were compliant with the European MLs of PCDD/Fs of 2 pg T fat and ALs of 1.75 pg TEQ/g fat, as well as with the European MLs of ndl-PCBs of 40 fat. Nonetheless, one sample from Lombardy exceeded the MLs of 4 pg TEQ/g fat lished for PCDD/Fs + dl-PCBs since contaminated by 6.28 pg TEQ/g fat, while the four samples from Lombardy were above the ALs of 2 pg TEQ/g fat established f PCBs ( Figure 2). When comparing cow milk contamination levels found in the present work with reported elsewhere in the world, a quite variable scenario emerged. For example, in Brazilian states, cow milk samples were found to have a PCDD/Fs contamination le 1.36 pg TEQ/g fat, which is approximately seven times higher than the concentration in this study [10]. The authors attributed the higher contamination levels to the incr proximity of farms to urban centers and the degree of industrialization in the area [1 contrast, cow milk samples collected from seven different regions of Chile showed co trations of PCDD/Fs + dl-PCBs up to 0.45 pg TEQ/g fat, two times lower than those When comparing cow milk contamination levels found in the present work with those reported elsewhere in the world, a quite variable scenario emerged. For example, in eight Brazilian states, cow milk samples were found to have a PCDD/Fs contamination level of 1.36 pg TEQ/g fat, which is approximately seven times higher than the concentration found in this study [10]. The authors attributed the higher contamination levels to the increased proximity of farms to urban centers and the degree of industrialization in the area [10]. In contrast, cow milk samples collected from seven different regions of Chile showed concentrations of PCDD/Fs + dl-PCBs up to 0.45 pg TEQ/g fat, two times lower than those found in this study [47]. The authors confirmed that higher contamination levels were observed in milk samples collected from more populated and industrialized areas [47].
When comparing data to those reported by Asian countries, it was observed that the PCDD/F levels of milk samples from traditional markets around Taiwan tended to be three-four times higher (0.59-0.89 pg TEQ/g fat) than those found in this study [48,49] but no significant differences in terms of dl-PCBs were identified [49]. Significantly lower PCDD/F + dl-PCB concentrations were recorded in milk samples from Malaysia (0.16 pg TEQ/g fat) [50], while concentrations of 0.13 ± 0.10 were reported in Chinese milk samples by Zhang et al. in 2015. However, since that value was expressed on a fresh weight rather than a lipid basis, it should be considered higher than those found in the present work [51].
In Europe, the EFSA reported a mean concentration of PCDD/Fs + dl-PCBs of 1.91 pg TEQ/g fat in milk and dairy products in 2012, which is higher compared to the contamination levels found in this study but includes contributions from milk-derived products (that usually have higher contamination levels than raw milk due to a larger amount of fat) [3]. In the same period, Marin et al. reported an average concentration of PCDD/Fs + dl-PCBs in milk and dairy products from Spain ranging from 0.25 pg TEQ/g fat to 2.13 pg TEQ/g fat [52].
Very high PCDD/F + dl-PCB contamination levels (up to 2-3 pg TEQ/g fat) were also reported in cow milk samples from Italy [7,53], as well as contamination levels up to 5.36 pg TEQ/g fat in buffalo milk from a contaminated site in southern Italy [54]. However, it is important to note that these data were calculated using the old toxic equivalency factors (TEFs) established in 1998 rather than those established in 2005 and currently employed [5]. Therefore, comparisons should be conducted with caution, especially since these surveys were conducted approximately 15 years ago and may not reflect the reduction trend observed in Europe in recent years.
Dioxin and PCB contamination levels in cow milk from Lombardy and Emilia-Romagna collected in this work (2018-2021) and the former monitoring plan (2012-2014) [37] were further compared. As can be observed from Table 2   The total TEQ values (pg TEQ/g fat) of both PCDD/Fs and dl-PCBs in the 2018-2021 milk samples from the two regions under study declined in direct proportion but less strongly than the actual mass values. Indeed, compared to the 2012-2014 outcomes, Lombardy and Emilia-Romagna saw PCB levels drop by 42% and 37%, respectively ( Table 2), suggesting that this drop was mostly due to less toxic PCDD/F and dl-PCB congeners (i.e., congeners having lower TEF values and, hence, a lower impact on total TEQ values).

Variation of Dioxins and PCB Congener Patterns
Since each dioxin and PCB congener bears the imprint of the contamination source (geographical location and industrial activities) and is characterized by a different toxicological effect, monitoring dioxin and PCB congener patterns in food is crucial to understand better their presence in the environment, contamination sources, and the potential risks to human health. Numerous studies have suggested that the chemical species of dioxins and PCBs may vary significantly depending on the foods being considered [3]. Regarding cow milk, dioxin, and PCB chemical species might differ based on a number of variables, including the feeding regimen, age, and breed of the animal, as well as the ambient conditions in which the cattle are grown, the sampling, and the analytical techniques employed for the quantification [55].
To provide complementary information, dioxin and PCB congener profiles in cow milk from the two Italian areas were examined for their variation based on the sample sites and investigated by plotting both the relative percentage contribution of each congener to the total PCDD/Fs + dl-PCBs TEQ values (Figure 3) and the absolute mass fraction concentrations (pg/g fat) (Figure 4).

Variation of Dioxins and PCB Congener Patterns
Since each dioxin and PCB congener bears the imprint of the contamination source (geographical location and industrial activities) and is characterized by a different toxicological effect, monitoring dioxin and PCB congener patterns in food is crucial to understand better their presence in the environment, contamination sources, and the potential risks to human health. Numerous studies have suggested that the chemical species of dioxins and PCBs may vary significantly depending on the foods being considered [3]. Regarding cow milk, dioxin, and PCB chemical species might differ based on a number of variables, including the feeding regimen, age, and breed of the animal, as well as the ambient conditions in which the cattle are grown, the sampling, and the analytical techniques employed for the quantification [55].
To provide complementary information, dioxin and PCB congener profiles in cow milk from the two Italian areas were examined for their variation based on the sample sites and investigated by plotting both the relative percentage contribution of each congener to the total PCDD/Fs + dl-PCBs TEQ values (Figure 3) and the absolute mass fraction concentrations (pg/g fat) (Figure 4).    As seen in Figure 3, PCB 126 was found to be responsible for 66% and 57% of the overall PCDD/Fs + dl-PCBs TEQ values in milk samples from Lombardy and Emilia-Romagna, respectively. These results contrast with those of Barone et al., according to whom cow milk showed a predominance of PCDD/Fs [56], but support earlier findings from our team, who found that PCB 126 was the predominant dl-PCB congener in milk samples in northern Italy [37,57]. Additionally, previous studies from other world regions identified PCB 126 and the other non-ortho PCBs as the congeners having the most impact on the total TEQ values in milk [47,[58][59][60][61].
The 2,3,4,7,8-PeCDF was found to be the chemical species with the second largest contribution (13% in samples from Lombardy and 9% in samples from Emilia-Romagna). As seen in Figure 3, PCB 126 was found to be responsible for 66% and 57% of the overall PCDD/Fs + dl-PCBs TEQ values in milk samples from Lombardy and Emilia-Romagna, respectively. These results contrast with those of Barone et al., according to whom cow milk showed a predominance of PCDD/Fs [56], but support earlier findings from our team, who found that PCB 126 was the predominant dl-PCB congener in milk samples in northern Italy [37,57]. Additionally, previous studies from other world regions identified PCB 126 and the other non-ortho PCBs as the congeners having the most impact on the total TEQ values in milk [47,[58][59][60][61].
The 2,3,4,7,8-PeCDF was found to be the chemical species with the second largest contribution (13% in samples from Lombardy and 9% in samples from Emilia-Romagna). On the other hand, two mono-ortho species, PCB 189 and PCB 105, showed a greater contribution to total TEQ values in milk samples from Emilia-Romagna (Figure 3).
PCDDs have typically been linked to agricultural chemicals such as chlorophenols, whereas more mixed PCDD/F profiles have been associated with both waste incineration and the usage of chlorinated agricultural pesticides [55]. Nonetheless, it can be challenging to accurately connect the PCDD and PCDF congener profiles of a specific food or feed to its sources because pollution sources might be numerous and varied.
As highlighted in Figure 3, milk samples from Lombardy also had a more pronounced pattern of two dl-PCB congeners, corresponding to PCB 105 and PCB 118 ( Figure 4B). In addition to particular PCB-containing technical formulations used in the past, pollution emissions from the steel sector, cement furnaces, and incineration facilities have been linked to the prevalence of these chemical species in food and feed [62]. Northern Italy hosted chemical industries producing chlorine compounds and actually houses a huge number of industrial facilities potentially responsible for the emission of different environmental pollutants [63]. Therefore, as already reported, the risk of contamination of local feed and food with PCDD/Fs and PCBs may be particularly high [37,57,64].
Differences in ndl-PCB patterns among the two regions were particularly evident for three congeners, namely PCB 138, PCB 153, and PCB 180 ( Figure 4C). Regardless of the sampling areas, these chemical species were the most abundant among the six ndl-PCB indicators, possibly because of their greater chlorination degree and longer persistence in the tissues of food-producing animals [55].

Screening-level Assessment via Deterministic Calculations
The UB exposure levels to PCDD/Fs + dl-PCBs due to the consumption of whole, semiskimmed, and skimmed cow's milk in all Italian age groups are presented in Table 3. The EWIs ranged from 0.02 to 15.23 pg TEQ/kg bw/week, with the lowest values calculated for mean elderly consumers ingesting skimmed milk and the highest values calculated for high infant consumers (95th percentile) ingesting whole milk. Globally, the EWIs of dioxins and PCBs were found to be higher across all population age groups due to the consumption of whole milk, while, regardless of the type of milk consumed, infants, toddlers, and children were found to be the most exposed population groups (Table 3).
A more thorough investigation revealed that for adolescents, adults, and the elderly, the EWIs from whole and semi-skimmed milk consumption were roughly the same (due to a compensation effect between the higher ingestions of semi-skimmed milk and the higher fat content of whole milk). On the other hand, despite the intermediate dietary consumption rates of skimmed milk, the EWIs associated with the consumption of this product were an order of magnitude lower ( Table 3).
The youngest segment of the European population was found to be much more likely to be exposed to dioxins and PCBs than the other groups by other authors. For example, Austrian children were reported to ingest 0.77 pg TEQ/kg bw/day of dioxins and PCBs through the whole diet, with milk and dairy contributing to 65% of the overall intake (EWI of 3.50 pg TEQ/kg bw/week) [65]. For French children, EWIs of 0.11-0.63 and 0.21-0.69 (0.77-4.41 and 1.47-4.83 pg TEQ/kg bw/week) were reported in relation to the consumption of cow milk and milk plus dairy products, respectively [66]. Table 3. Estimated weekly intake (EWI, pg TEQ/kg bw/week) to PCDD/Fs + dl-PCBs of different Italian population age groups due to the chronic mean and high (95th percentile, in square brackets) consumption of whole (3.5 g fat/100 g), semi-skimmed (1.8 g fat/100 g) and skimmed (0.3 g fat/100 g) cow milk containing average concentrations of 0.78 pg TEQ/g fat of PCDD/Fs + dl-PCBs.

Population Group
Type of Milk As for Italian children, the whole diet was estimated to be responsible for a weekly intake of 1.98-4.98 pg TEQ/kg bw/day, with milk alone contributing to 10% (EWI of 1.38-3.48 pg TEQ/kg bw/week) and, when combined with dairy products, to 34% (EWI of 4.71-11.85 pg TEQ/kg bw/week) of the weekly ingestions of these contaminants [67]. Similarly, in another study, 1.50-2.64 pg/TEQ/kg bw/week of dioxins and PCBs were estimated to be taken by Italian children through the consumption of milk and dairy products [56]. Hence, based on the calculated EWIs and the conventional deterministic risk characterization analysis summarized in Figure 5, a potential health risk concern related to the consumption of whole and semi-skimmed milk emerged, with infants, toddlers, and children being the populations at higher risk. Indeed, average intakes significantly exceeded the toxicological reference values, being 125-320% (95th percentile: 390-762%) and 85-175% (95th percentile: 192-409%) of the TWI for infants and toddlers, respectively ( Figure 5). As for children, the contribution of whole and semi-skimmed milk consumption to the TWI was lower (41-82%; 95th percentile: 95-213%) but still concerning, given that many foods, besides milk, provide significant amounts of dioxins and PCBs ( Figure 5).
Our risk characterization provides more alarming results compared to those previously reported in Italy, according to which cow milk contributed roughly 20-30% of the TWI [37,53]. Nevertheless, it should be noted that a direct comparison can be misleading since these studies used old dioxin TEF values or compared the exposure level to the Tolerable Daily Intake (TDI) of 2 pg WHO-TEQ/kg bw, which is seven times lower than the actual TWI. Additionally, the conservative nature and simplifications of the deterministic exposure assessment, which means lacking information about the variability and uncertainty in both food consumption statistics and contamination values, typically lead to incomplete or inaccurate (i.e., overestimated) exposure levels, giving information on the possibility but not the probability of risk.
sumption of whole and semi-skimmed milk emerged, with infants, toddlers, and children being the populations at higher risk. Indeed, average intakes significantly exceeded the toxicological reference values, being 125-320% (95th percentile: 390-762%) and 85-175% (95th percentile: 192-409%) of the TWI for infants and toddlers, respectively ( Figure 5). As for children, the contribution of whole and semi-skimmed milk consumption to the TWI was lower (41-82%; 95th percentile: 95-213%) but still concerning, given that many foods, besides milk, provide significant amounts of dioxins and PCBs ( Figure 5). Our risk characterization provides more alarming results compared to those previously reported in Italy, according to which cow milk contributed roughly 20-30% of the TWI [37,53]. Nevertheless, it should be noted that a direct comparison can be misleading since these studies used old dioxin TEF values or compared the exposure level to the Tolerable Daily Intake (TDI) of 2 pg WHO-TEQ/kg bw, which is seven times lower than the actual TWI. Additionally, the conservative nature and simplifications of the deterministic exposure assessment, which means lacking information about the variability and uncertainty in both food consumption statistics and contamination values, typically lead to incomplete or inaccurate (i.e., overestimated) exposure levels, giving information on the possibility but not the probability of risk.

Refined Assessment via Probabilistic Modelling
Monte Carlo simulations were run to estimate in a more accurate and realistic way the distribution of EWIs and the contribution of milk consumption to the TWI and, therefore, to quantify the probability of the health risk occurring. Only those exposure scenarios of particular concern resulting from the deterministic assessment were simulated. Moreover, since the imputation method for the treatment of data below the LOQ might have a large effect on the forecasted outcomes, simulations were run three separate times using UB, MB, or LB concentrations one at a time. Figure 6 presents the variability in the form of distribution frequencies of EWIs at the MB scenario for the most vulnerable population consuming semi-skimmed and whole milk. Average EWIs were estimated to range from 5.74 ± 7.77 pg/TEQ/kg bw/week for infants consuming whole milk to 0.74 ± 0.97 pg/TEQ/kg bw/week in children consuming semi-skewed milk. The median (P50) MB values were significantly lower when compared to the mean values ( Figure 6) because of the right-skewness of the distribution frequency curve, which resulted as such because of the right-skewness of both consumption data

Refined Assessment via Probabilistic Modelling
Monte Carlo simulations were run to estimate in a more accurate and realistic way the distribution of EWIs and the contribution of milk consumption to the TWI and, therefore, to quantify the probability of the health risk occurring. Only those exposure scenarios of particular concern resulting from the deterministic assessment were simulated. Moreover, since the imputation method for the treatment of data below the LOQ might have a large effect on the forecasted outcomes, simulations were run three separate times using UB, MB, or LB concentrations one at a time. Figure 6 presents the variability in the form of distribution frequencies of EWIs at the MB scenario for the most vulnerable population consuming semi-skimmed and whole milk. Average EWIs were estimated to range from 5.74 ± 7.77 pg/TEQ/kg bw/week for infants consuming whole milk to 0.74 ± 0.97 pg/TEQ/kg bw/week in children consuming semi-skewed milk. The median (P50) MB values were significantly lower when compared to the mean values ( Figure 6) because of the right-skewness of the distribution frequency curve, which resulted as such because of the right-skewness of both consumption data and contamination levels used as input values for simulations. This condition led to greater robustness of EWIs at the left tail of the curve (lower percentiles). On the contrary, a greater degree of uncertainty is associated with the interpretation of EWIs beyond the 90th percentile since they are not supported by an adequate number of observations. In addition, from the comparison of EWIs at different 10% increment percentiles (Table 4), it was observed that there was a significant amount of uncertainty due to the method for the imputation of data below the LOQ. Indeed, the MB estimates were roughly 15% higher or lower than the LB or UB estimates, respectively. All these findings suggest a possible overestimation of the higher percentile intakes.
Similar considerations also apply to the outcomes resulting from the Monte Carlosimulated risk characterization associated with chronic exposure to dioxins and PCBs. The MB 95th percentile contribution percentages of milk ingestion to the TWI ranged, on average, from a minimum of 22% for children ingesting semi-skimmed milk to a maximum of 169% for infants ingesting whole milk (Table 5), being significantly lower than those resulting from the deterministically calculated worst case scenario ( Figure 5). Nevertheless, as it can be observed from Table 5, the consumption of whole milk by infants was already above the TWI at the MB 40th percentile, indicating that this population group is at particular risk due to the intake of dioxins and PCBs. At the 80th percentile, only the consumption of semi-skimmed milk by children was below the TWI. a greater degree of uncertainty is associated with the interpretation of EWIs beyond the 90th percentile since they are not supported by an adequate number of observations. In addition, from the comparison of EWIs at different 10% increment percentiles (Table 4), it was observed that there was a significant amount of uncertainty due to the method for the imputation of data below the LOQ. Indeed, the MB estimates were roughly 15% higher or lower than the LB or UB estimates, respectively. All these findings suggest a possible overestimation of the higher percentile intakes.     To quantify the uncertainty around when each exposure scenario can occur, the cumulative probability distributions of milk intake percentage contributions to the TWI for infants, toddlers, and children were plotted in Figure 7. Ultimately, the probability for infants to exceed the TWI was 76% for the ingestion of whole milk and 35% for the ingestion of semi-skimmed milk, which means that three out of four and one out of three infants are at a very high risk of chronic toxicity. Moreover, 56% and 26% of toddlers, as well as 22% and 8% of children consuming whole and semi-skimmed milk, respectively, can be exposed to dioxin and PCB levels above the toxicological reference values (Figure 7). Such high percentages were already reported for infants and toddlers through the whole diet [68], but our findings might have more severe implications from a health perspective, given that they are linked to the sole consumption of cow milk, while the diet includes many other sources of these contaminants such as meat, fish, and eggs [3].
Despite these results and the high confidence in probabilistic predictions, there are a few limitations related to the applied methodology that need to be mentioned. In particular, regardless of the substitution method applied, the precision of the measured contaminants may have been biased by the left censorship of these data. Dietary intakes of milk by the general Italian population may not accurately reflect those of the inhabitants of northern Italy due to potential regional variations in dietary habits. Moreover, it was assumed that the entire amount of cow milk consumed originated from local farms in Lombardy and Emilia-Romagna, while commercial products can usually be a mixture of milk from different European (or non-European) countries. The unavailability of data concerning the distribution of the population's body weight might have been another important limiting factor, which, during probabilistic calculations, has not made it possible to evaluate its impact on the output of the exposure assessment. Finally, the potential variation in the bioavailability of dioxins and PCBs due to milk processing and/or consumer age was not taken into consideration.
22% and 8% of children consuming whole and semi-skimmed milk, respectively, can be exposed to dioxin and PCB levels above the toxicological reference values (Figure 7). Such high percentages were already reported for infants and toddlers through the whole diet [68], but our findings might have more severe implications from a health perspective, given that they are linked to the sole consumption of cow milk, while the diet includes many other sources of these contaminants such as meat, fish, and eggs [3].

Conclusions
Despite confirming a significant downward trend in dioxin and PCB concentrations within European food and feed chains, this study raises important questions about the potential health risks associated with dietary long-term exposure to these contaminants for the youngest population. By applying probabilistic techniques, a more accurate and realistic picture of exposure levels to dioxins and PCBs and associated risks were achieved, which allowed estimating that the sole consumption of whole cow milk (albeit broadly compliant with maximum EU limits) may lead to exceeding the tolerable weekly intake in approximately 76, 56, and 22% of infants, toddlers, and children, respectively. This output, which is mainly the consequence of the recent seven-fold reduction in the tolerable weekly intake of dioxins and PCB, clearly worsens the overall risk profile of the young population and underscores the need for a multi-pronged approach involving regulatory measures and public education to reduce dietary exposure further. In this setting, monitoring and surveillance systems represent an essential tool to identify contamination sources and track human exposure. On the other hand, effective communication of both risks and benefits associated with the consumption of milk by vulnerable consumers is necessary. This may include promoting a varied diet and alternating with infant formulae, follow-up formulae, and other baby foods for which stricter limits of dioxins and PCBs are in force.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/foods12091869/s1, Table S1: Limit of quantification (LOQ) and recovery values of the 35 measured PCDD/F and PCB congeners. Table S2: Mean concentrations of the 35 measured PCDD/F and PCB congeners in milk samples from the two Italian regions. Table  S3 Funding: This research did not receive any specific funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The dataset generated during the current study will be made available upon reasonable request.