Chemometric Optimization of Solid-Phase Extraction Followed by Liquid Chromatography-Tandem Mass Spectrometry and Probabilistic Risk Assessment of Ultraviolet Filters in an Urban Recreational Lake

Solid-phase extraction (SPE) of eleven ultraviolet filters (UVFs): benzophenone-1 (BP-1); benzophenone-3 (BP-3); benzophenone-4 (BP-4); isoamyl p-methoxycinnamate (IAMC), homosalate (HMS); 4-hydroxybenzophenone (4-HB); 4-methylbenzylidene camphor (4-MBC); octocrylene (OC); octyl dimethyl-p-aminobenzoate (OD-PABA); 2-ethylhexyl-4-methoxycinnamate (EHMC); and avobenzone (AVO), has been optimized using Plackett-Burman design, Box-Behnken design, and Derrindzer desirability function. Of the six SPE variables studied, the most influencing is the type of eluent followed by pH and the methanol content in the rinsing solvent. A method with good analytical performance was obtained by applying optimal SPE conditions and liquid chromatography-tandem mass spectrometry (LC-MS/MS), with the method detection limit ranging from 0.1 to 5 ng/L, recovery from 44% to 99%, and relative standard deviation (RSD) within 19%. This method was used to analyze the content of UVFs in an urban lake (Sava Lake, Serbia). UVFs occurrence, geostatistical distribution, and associated environmental risk are highly dependent on recreational activities. The average concentrations of UVFs ranged from 0.3 to 113 ng/L, and the most present substance was EHMC, followed by 4-MBC and BP-3. The spatial distribution of the risk quotient (RQ = 0.04–1.7) inside the lake is highly correlated with the number of people bathing and swimming. Human exposure through the dermal pathway is higher than ingestion for most UVFs. Monte Carlo simulation of probabilistic risk assessment estimated the percentile P10, P50, P90 of 12.7; 17.3; 47.5 and 20.1; 27.6; 77.5 ng/kg∙day for total human exposure of adults and children, respectively. The sensitivity analysis revealed that the health risk estimate depends mostly on the content of EHMC, HMS, and 4-MBC, while the most influential exposure variables were human body weight and skin surface area. There is no serious concern to human health due to UVFs in the short term; however, a high ecological risk in some parts of the lake is estimated.


INTRODUCTION
Ultraviolet filters (UVF) are used to protect the human skin from harmful UV radiation and are added to various skin treatment products (Kunze et al., 2021;Kwon and Choi, 2021). However, some UVFs, especially if applied in large quantities, can cause side effects, such as various allergies, dermatitis, and adverse effects on hormonal receptors, nervous system, and reproductive functions (Huang et al., 2021).
Concentrations of individual UVFs in sunscreen preparations are at the level of several percent, and the European Commission (EC No 1223) has set limit/reference values for individual organic UVFs in the range of 2%-10% (except for drometrizole trisiloxane 15%). UVFs are used daily, and their use is exceptionally high during the summer months (Fagervold et al., 2019), which leads to the release of vast amounts of these substances into the environment (Ma et al., 2016).
Unfortunately, UVFs have been proved to pose a high risk for many aquatic organisms (Carve et al., 2021). Different UVFs affect the environment differently (Labille et al., 2020). Some UVFs tend to accumulate in living cells, affecting reproduction, and increasing mortality in benthic organisms, small fish, and sizeable newborn fish (Carve et al., 2021). This suggests that the presence of UVFs in the water should be considered with great concern.
A significant proportion of UVFs enters the environment through sewage effluents (Wang et al., 2021), which may further contaminate the surface water and groundwater (Li et al., 2020a). In addition, one of the ways of polluting the environment water with UVFs is rinsing the skin during swimming (Teo et al., 2015). Therefore, the presence of UVF needs to be investigated in water at the wastewater treatment plants and in the vicinity of beaches.
Sava Lake in Belgrade (Serbia) is a typical urban lake that is used for recreation such as swimming, sailing, fishing, and diving. In the summer months, there are several thousand people at the lake daily, and it is expected that a large amount of UVFs will be released into the lake water. However, the recreational activities are not evenly distributed along the lake, i.e., there are parts of the lake where there is a big crowd with people sunbathing and swimming and parts of the lake without such activities. Consequently, the different spatial distribution of potential water contaminants is expected.
Risk assessment studies have shown an increased ecological risk due to UVFs in recreational locations and that it depends on the intensity of application of sunscreen products (Teo et al., 2015). Keeping in mind that the concentrations of different UVFs in waters are different and each UVF has a different distribution, the probabilistic approach to risk assessment in swimming water has recently become irreplaceable (Astuti et al., 2022;Dehghani et al., 2022). Compared to the deterministic one, this approach gives a complete picture of the risk estimate, considering that individual UVFs concentrations differ temporarily and spatially (Tsui et al., 2015). In addition to the distribution of UVF concentrations, the distribution of the exposure parameters influencing the risks should not be neglected.
Several analytical methods can be used to quantify the UVFs traces in water. When UVF concentrations are very low, the most optimal technique for sample preparation is solid-phase extraction (SPE) (Allinson et al., 2018;Cadena-Aizaga et al., 2020). This technique enables high degrees of preconcentration, which is why it is the most commonly used technique for analyzing trace level organics in water. In order to get the best analytical performance, the SPE process needs to be optimized (Silva et al., 2013;Chávez-Moreno et al., 2018). As for the instrumental measurement itself, two techniques are predominantly used to measure UVFs. One is gas chromatography with tandem mass spectrometry (GC-MS/MS) (da Silva et al., 2015;Huang et al., 2016;Ramos et al., 2019), which has a huge potential for identification and selectivity in separation. Another technique is liquid chromatography with tandem mass spectrometry (LC-MS/ MS) (Li et al., 2020b;Chiriac et al., 2020;Cadena-Aizaga et al., 2022), which is characterized by exceptional sensitivity and is very suitable for the analysis of lowvolatile and degradable UVFs.
The main objectives of this work were to optimize the SPE sample preparation prior to LC-MS/MS, investigate the occurrence and spatial distribution of different UVFs in recreational water of an urban city lake, and assess potential environmental risk using a probabilistic approach.

Study Area and Sampling
The Sava Lake (Figure 1), better known as Ada Ciganlija, is a tributary of the Sava River. The lake was formed about 50 years ago. The length of the lake is about 4.4 km, the average width is 250 m, and the average depth is 4.5 m (Trbojević et al., 2021). The lake volume is about 4 million m 3 , and the area is approximately 1.0 km 2 . On the left side, there is a small fenced pool connected to the Sava River, from which water is pumped into the Sava Lake, while on the right side, water is discharged from the lake into the Čukarica stream. This artificially ensures water flows through the lake.
Sampling was performed in a one-time sampling campaign conducted in August 2021, at the peak of the recreational season. The population of swimmers/bathers on Sava Lake at the time of sampling was about a thousand people. A total of 28 grab samples were taken along the coast and in the middle of Sava Lake. The samples were taken 10 cm below the water surface, stored at 4°C, and analyzed within 24 h. All water samples were collected in amber glass bottles (volume 0.5 L).

Sample Preparation
A two-step statistical approach was employed to optimize several variables that affect the SPE sample preparation of UVFs. The variables, such as the type of eluent, percentage of methanol in the rinsing solution (%MeOH), initial pH value, the sample volume (V s ), the rinsing solution volume (V r ), and eluent volume (V e ), were screened by applying Plackett-Burman design (PBD) and were optimized further by central composite design (CCD). The variable and their coded values for PBD are listed in Table 1. PBD experiments, consisting of 15 runs, have been performed according to Supplementary Tables S1, S2, while the subsequent CCD experimental runs are listed in Supplementary Table S3.
The effect of pH on the UVFs extraction from sample solution was studied within the range of 2-8, using 0.1 M HCl and 0.1 M NaOH for pH adjustment. Spiked samples were made by dissolving standard mixture solutions of 1.0 μg/ml each UVF in water. After adjusting the initial pH values of the samples, the conditioning of the SPE cartridge was started. A 5 ml volume of the eluent was applied to the cartridges, followed by 5 ml of deionized water whose pH was adjusted to the same value as the pH of the sample to be applied. The sample application was made at approximately 1.0 ml/min flow rate. After applying the sample, the cartridges were washed with a rinsing solution and dried for 10 min. The elution was performed with the selected type of the eluent and its optimized volume. The eluted extracts were evaporated under a stream of nitrogen and reconstituted with 0.5 ml of methanol. All samples were filtered with a 0.22 μm PTFE syringe filter.

Liquid Chromatography-Tandem Mass Spectrometry Measurement
A Thermo Fisher Scientific (Waltham, United States) LC-MS/MS system comprised of a Dionex Ultimate 3000 liquid chromatograph fitted with an Accucore ™ aQ C-18 column (10 cm length, 2.1 mm I.D., particle size 2.6 µm) and LTQ XL linear ion trap mass spectrometer operating in the electrospray ionization (ESI) mode. Mass spectra of analytes were processed by an Xcalibur Ver. 3.0 software. A 10 µL volume of the SPE prepared extract was injected into the LC-MS/MS system. LC mobile phase composition and MS/MS operating parameters are given in Supplementary Tables S4, S5.

Data Analysis
First, the experimental results were processed by descriptive statistics, Pearson correlation, and the Ryan-Joiner normality test. Then, the log-transformed concentration data were subjected to hierarchical cluster analysis (HCA) using Ward linkage and Euclidean distance. Finally, geostatistical distribution analysis was based on the ordinary kriging  method, and results were mapped in the polygon representing Sava Lake. The ecological risk associated with UVFs in the Sava Lake was assessed based on the risk quotient (RQ), which is the ratio between a measured concentration (C) and a predicted no-effect concentration (PNEC), as described in Eq. 1: Human exposure to UVFs through ingestion and dermal pathways has been evaluated by following the US EPA SWIMODEL exposure assessment model (Li et al., 2020b;Dehghani et al., 2022). Accordingly, the calculation formulas are given in Eqs 2, 3, while the exposure parameters and their values used for risk assessment are listed in Table 2 and  Supplementary Tables S6, S7.
Monte Carlo simulation and sensitivity analysis have been made with the following run preferences: the number of trials is 1000, the confidence level is 95%, and the initial seed value of 999.

Plackett-Burman Design Screening
Since all UVFs are simultaneously subjected to the SPE process, additional efforts should be made to obtain an optimal response that would include all individual responses. Thus, considering that we have 11 UVFs, an aggregate response that includes the individual responses for all analytes simultaneously is used here.
The most commonly used cumulative response in analytical measurements is the desirability function (D) proposed by Derringer and Suich (1980), which represents the geometric mean estimated as a product of the individual desirabilities (d i ) weighted by their individual coefficients (r i ) (Eq. 4). In this case, no difference in weights for individual responses was assumed, so the D values were calculated taking into account all UVFs equally. Therefore, the values for d i have been maximized and ranged from 0 to 1, while the r i values were fixed (r i = 1).
The results of the screening design are given in the form of a Pareto plot in Figure 2. It shows that the parameter of the most significant influence is the type of SPE eluting solvent. Besides, a considerable influence is also noticed for %MeOH and pH, while other parameters have much less impact on the SPE process. The main effects plot for the tested analytes (Supplementary Figure S1) revealed that in addition to the type of solvent, recoveries are positively affected by increasing the parameters V e and V s . Conversely, the increase in the values of the parameters %MeOH, pH, and V r affect the response negatively.
In the second step, ethyl acetate (EtOAc) was selected as an eluent, and the percentage of methanol (%MeOH) in the rinsing solvent and the sample pH value were selected as variables for the CCD optimization.

Central Composite Design Optimization
The initial pH of the sample and %MeOH have been optimized according to CCD in the experimental domain for pH = 2-8 and %MeOH = 2%-8%. Figure 3 shows the response surface plot for optimization of pH and %MeOH. The response was fitted to a second-order polynomial model represented by the following regression equation: Final SPE conditions were obtained by optimizing the D value ( Figure 3). As a result, a maximum D value of 0.67 was obtained for %MeOH of 4.8% and pH = 5.1. Finally, the global optimum for SPE of all UVFs was selected for the following values of variables: the extraction volume of 150 ml of aqueous sample, pH of the sample set at 5.1, and 15 ml aqueous solution with 4.8%MeOH to rinse the cartridge, and 15 ml EtOAc for elution. These conditions were used to prepare water samples from Sava Lake for the analysis of UVFs by the LC-MS/MS method.

Validation of SPE-LC-MS/MS Method
The optimized SPE sample preparation was accompanied by LC-MS/MS measurement. The LC-MS/MS method was previously optimized, which included optimizing the HPLC mobile phase (Vasiljević et al., 2004), tuning the MS instrument, finding optimal m/z and values, and CID value optimization (Grujić et al., 2009). The optimal LC-MS/MS parameters are given in Supplementary Tables S4, S5. Typical LC-MS/MS In order to assess the analytical performances of the method, the method detection limits (MDL), average recovery (R), and relative standard deviation (RSD) were determined. Clean surface water, in which the analyzed UVFs were not detected, was used to make spiked samples and construct matrix-matched calibration lines at nine concentration levels for the UVFs mixture. The calibration range was 0.1-300 ng/L. For each UVF, the linear correlation coefficient was greater than 0.995. The obtained MDL, R, and RSD values are listed in Supplementary Table S8. MDLs, estimated as three times the signal-to-noise ratio, range from 0.01 to 4.4 ng/L. Recovery studies were performed at 1.0, 10, and 100 ng/L concentration levels. Average recoveries for 11 analytes range from 44% to 99%. Six replicates were used to investigate repeatability at the mentioned concentration levels, and RSDs values within ±19% were obtained.
Pearson correlation matrix (Supplementary Table S9) reveals a strong association between UVFs. Only OD-PABA concentrations are not related to other UVFs. The highest correlation coefficient of 0.969 was found for the pair EHMC and AVO, while the smallest one was estimated for the pair BP-1 and OD-PABA (0.102). In general, it can be concluded that there is a large intercorrelation between UVFs. Figure 5 shows a dendrogram obtained by hierarchical cluster analysis (HCA). The relationship between the groups of tested samples and UVFs is obvious. Four clusters (I, II, III, and IV) were identified. This grouping can be attributed to the different spatial distribution of sampling points and the difference in the amount of UVFs in sunscreen products and how they are  applied. Namely, sunscreen products are usually composed of a combination of UVFs to ensure the broadest protection.

Deterministic Approach
RQs for individual UVFs were obtained by dividing C i obtained in this study by PNECs values (Supplementary  Table S10), which are collected after the literature survey (Lukić et al., 2021). The risk classification was based on risk ranking criteria in which RQ < 0.01; 0.01 < RQ < 0.1; 0.1 < RQ < 1; and RQ > 1, are "Unlikely to pose risk"; "Low risk"; "Medium risk" and "High risk", respectively (Ma et al., 2016).
Contributions of individual UVFs to total RQ are presented in Figure 6. It is noticeable that the greatest contribution to RQ is made by the presence of EHMC, followed by 4-MBC and BP-1. The share of other UVFs in total RQ can be almost neglected. The estimated minimum, average and maximum RQs values are 0.4, 0.35 and 1.7, respectively. There are two groups of samples in which the RQ values exceed 1. The first group comprises samples that highly exceed one and are marked 13-15, while the second group includes samples numbered 12, 27, and 28, slightly above 1.
In conclusion, both sample groups represent locations with a high ecological risk.
The geostatistical distribution of risk in Sava Lake was made based on calculated RQs from the measured values, which were used to perform ordinary kriging interpolation on the entire lake area (Radomirović et al., 2020). Figure 7 shows that in several places close to the lake beach, the RQs are elevated in comparison to the middle of the lake. Several hotspots have been identified (A, B, C, D), which correlate with the intensity of human activities at the site.
In area A, samples 12, 13, 14, and 15 are located, while samples 27 and 28 are in area B. These samples show RQ ˃ 1 ( Figure 6) and are classified in clusters III and IV in the dendrogram ( Figure 5). Samples 19 and 20 are in area D, and sample 22 represents sampling point C. This distribution reflects the intensity of human activity at the lake. Most people and swimmers/bathers are located on the beach and in the water in area A, then an area at point B. At locations C and D, there are beaches with significantly fewer people. In the middle of the lake, as well as at the left and right edges, water contamination with UVFs is significantly lower.

Monte Carlo Simulation and Sensitivity Analysis
Since UVFs concentrations are not uniformly distributed across the lake, a probabilistic approach to human exposure is expected to provide a more realistic assessment than a point estimate (Astuti et al., 2022;Dehghani et al., 2022). Given that the exposure parameters also have their distributions, individual risks and cumulative risk will also follow some distributions. Namely, the exposure parameters AT and BI are fixed points, SF and ST have triangle distribution, while BW and SA follow a log-normal distribution. PNEC, Log Kow, and MW are analyte-specific fixed-values ( Table 2,  Supplementary Tables S6, S7). Concentrations of individual UVFs have been fitted to estimate the best distribution. All UVFs concentrations follow a log-normal distribution, except BP-4, which is best fitted with the Webuill curve.
Monte Carlo simulation of human exposure was applied to UVFs in Sava Lake. According to the SWIMMER model, humans may be exposed to UVFs through inhalation, digestion, and dermal pathway. While inhalation is less likely to be a significant route due to the low volatility of UVFs (Lu et al., 2017), the other two are expected to be the important exposure routes of UVFs to humans. The accidental ingestion and dermal absorption of UVFs were calculated from Eqs 2, 3. The exposure parameters listed in Table 2 were used to estimate the human dose for adults and children (age   Supplementary Table S11. Human exposure through the dermal pathway poses a higher risk compared to ingestion for most UVFs, except for BP-1, BP-3, BP-4 and 4-HB. Figure 8 shows the human exposure distribution for adults and children and the corresponding sensitivity risk analysis. The histograms for the total exposure dose show the simulated ΣD a and ΣD c in 90% (the 90 of the estimated percentile) for the adult and children groups of 47.5 and 77.5 ng/kg·day, respectively. This indicates that children are much more at risk due to UVFs exposure than adults because of their lower body-weight/skin surface ratio. Anyway, the lack of reference dose (R f D) values prevents the explicit expression of the health risk index. Instead, the UVFs exposure can be compared to similar substances, i.e., parabens for which there is R f D of 10 mg/kg·BW/day (Li et al., 2021) proposed by the Joint Expert Committee on Food Additives (JECFA) of Food and Agriculture Organization of the United Nations (FAO) and World Organization (WHO). Based on this comparison, it can be concluded that there is a negligible risk to human health for the estimated UVFs doses (ng/kg·day) (Li et al., 2020B). Thus, the concentrations of these compounds in Sava Lake water can be considered acceptable, except at a few hotspot locations. Anyhow, there are concerns that UVFs will affect human health in the long term either through their influence or from harmful organic substances that may arise as products of their degradation.

CONCLUSION
SPE may be efficiently step-wise optimized using screening and response surface designs. It is a reliable sample preparation technique, which, coupled with LC-MS/MS, represents a suitable analytical tool for UVFs in surface water. Among 11 popular UVFs, EHMC, 4-MBC, and BP-3 were most abundant in   S, skewness (unitless); K, kurtosis (unitless).
FIGURE 8 | Probability and sensitivity analysis of different groups exposed to UVFs: (A) adults; (B) children.
Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 916916 the studied recreational lake. It was shown that recreational activities greatly influenced the occurrence of UVFs, their spatial distribution, and potential environmental risk. UVFs are significantly present near the coast, i.e., close to the beach. Nevertheless, the RQ values exceeded one at several locations in the lake. The incidental ingestion is a minor route to human exposure to most UVFs. Monte Carlo simulation of risk estimated 63% higher percentiles for total human exposure of children compared to adults. In addition, the sensitivity identified the most influential variables as the content of EHMC, HMS, and 4-MBC, along with two exposure parameters, human body weight and skin surface area. In general, except for the high ecological risk in some parts of the lake, UVFs in studied recreational lake water posed no concern to human health in the short term.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JL: methodology, investigation, and writing-original draft. JR: resources and validation. ML: formal analysis and visualization. TĐ: supervision, project administration, and funding acquisition. AO: conceptualization, software, writing-review and editing.