Using Design of Experiments to Optimize a Screening Analytical Methodology Based on Solid-Phase Microextraction/Gas Chromatography for the Determination of Volatile Methylsiloxanes in Water

Volatile methylsiloxanes (VMSs) constitute a group of compounds used in a great variety of products, particularly personal care products. Due to their massive use, they are continually discharged into wastewater treatment plants and are increasingly being detected in wastewater and in the environment at low concentrations. The aim of this work was to develop and validate a fast and reliable methodology to screen seven VMSs in water samples, by headspace solid-phase microextraction (HS-SPME) followed by gas chromatography with flame ionization detection (GC-FID). The influence of several factors affecting the extraction efficiency was investigated using a design of experiments approach. The main factors were selected (fiber type, sample volume, ionic strength, extraction and desorption time, extraction and desorption temperature) and optimized, employing a central composite design. The optimal conditions were: 65 µm PDMS/Divinylbenzene fiber, 10 mL sample, 19.5% NaCl, 39 min extraction time, 10 min desorption time, and 33 °C and 240 °C as extraction and desorption temperature, respectively. The methodology was successfully validated, showing low detection limits (up to 24 ng/L), good precision (relative standard deviations below 15%), and accuracy ranging from 62% to 104% in wastewater, tap, and river water samples.


Introduction
Siloxanes are organic compounds with a backbone of alternating atoms of silicon (Si) and oxygen (O), with organic groups (e.g., methyl, ethyl, phenyl, etc.) attached to the silicon atoms [1]. The most significant representatives are the low molecular volatile methylsiloxanes (VMSs), which are commonly incorporated in daily products such as personal care products (PCPs). Considering their chemical structure, they may be classified as linear (lVMSs) or cyclic (cVMSs) compounds [2]. Linear and cyclic methylsiloxanes are usually expressed as Ln and Dn, with "n" representing the number of silicon atoms.
Due to their massive use, VMSs are continuously discharged down-the-drain, reaching wastewater treatment plants (WWTPs). Conventional WWTPs are not prepared to degrade them. Thus, VMSs may accumulate on sewage sludge or be discharged into the environment (watercourses) via final effluent [3,4]. Cyclic siloxanes (specially D4, D5, D6) are usually detected at higher levels than linear, possibly due to their higher production 2.1.1. Screening Design Several parameters may influence the HS-SPME extraction. Therefore, a screening design considering seven main factors (Xi) was implemented. Two levels (defined as −1 and +1 for the lower and upper limit, respectively) were set for each of the seven main factors, as presented in Table 1: X1-ionic strength (0-20%, w/v), X2-extraction time (5-45 min), X3-desorption time (1-10 min), X4-extraction temperature (25-80 • C), X5desorption temperature (200-250 • C), X6-type of fibers (100 µm polydimethylsiloxane (PDMS), 65 µm PDMS/Divinylbenzene (DVB)) and X7-sample volume (5-10 mL). A 2 7-3 fractional factorial screening design was chosen (resolution IV), leading to a total of 16 experiments, which corresponds to one eighth of the full factorial design with 128 experiments ( Table 2). All the extractions were tested using a concentration of 1 µg/L for both VMSs and internal standard (M4Q). The main results are presented in Figure 1. The model was applied to all compounds and its suitability determined. All responses were adjusted to a quadratic model with a R 2 > 0.94. The main effects were determined by the F-probability, where F-probability ≤ 0.05 represents a significant effect and if 0.05 < F-probability ≤ 0.1, a relative effect may be considered. Since the main goal was to maximize the response areas for all compounds, a desirability function was applied in order to obtain the best compromise value for all variables. represents a significant effect and if 0.05 < F-probability ≤ 0.1, a relative effect may be considered. Since the main goal was to maximize the response areas for all compounds, a desirability function was applied in order to obtain the best compromise value for all variables. As presented in Figure 1, all main factors were significant for at least one compound and the extraction temperature (X4), fiber type (X6), and sample volume (X7) were the factors with stronger effects on the response (F-probability < 0.0001). The fiber type is a discrete variable, and therefore to proceed with the model, it had to be defined in advance. A PDMS/DVB fiber was selected since it provided higher responses for all compounds. In this type of fiber, the extraction mechanism is based on adsorption into pores, which is a stronger and more efficient mechanism than absorption, making this fiber more suitable for trace analysis. Analyzing the screening design results, the larger the sample volume, the greater the instrumental response. Therefore, the maximum sample volume tested, 10 mL, was chosen (half of the vial capacity). These results are consistent with two other studies found in the literature that use HS-SPME for the analysis of VMSs in water, although GC-MS was used in both [12,13].

Central Composite Design
After selecting the most important factors using a screening design, those selected parameters were optimized by a central composite design (CCD). Considering that five out of the seven factors were relevant (ionic strength, extraction and desorption time, and extraction and desorption temperature), a total of 26 experimental runs were performed, including six assays in the center of the cubic domain (pattern 00000), resulting in a total of 32 experiments. The CCD factor values and experiments are summarized in Supplementary Tables S1 and S2. The experimental data were fitted to a second-order polynomial equation and the model coefficients were determined based on a least-square regression analysis. The equations and the model suitability were assessed using the ANOVA approach and the main results are presented in Supplementary Table S3 (only the factors that significantly affect the response are shown). The R 2 values ranged from 0.80 to 0.97, which indicates a reasonable relationship between the experimental data and the fitted model. Models for most compounds show a F-probability below 0.1, indicating that variations in the response are associated with the models, rather than with the experimental error. However, M4Q and L5 show a higher F-probability, indicating a strong influence As presented in Figure 1, all main factors were significant for at least one compound and the extraction temperature (X4), fiber type (X6), and sample volume (X7) were the factors with stronger effects on the response (F-probability < 0.0001). The fiber type is a discrete variable, and therefore to proceed with the model, it had to be defined in advance. A PDMS/DVB fiber was selected since it provided higher responses for all compounds. In this type of fiber, the extraction mechanism is based on adsorption into pores, which is a stronger and more efficient mechanism than absorption, making this fiber more suitable for trace analysis. Analyzing the screening design results, the larger the sample volume, the greater the instrumental response. Therefore, the maximum sample volume tested, 10 mL, was chosen (half of the vial capacity). These results are consistent with two other studies found in the literature that use HS-SPME for the analysis of VMSs in water, although GC-MS was used in both [12,13].

Central Composite Design
After selecting the most important factors using a screening design, those selected parameters were optimized by a central composite design (CCD). Considering that five out of the seven factors were relevant (ionic strength, extraction and desorption time, and extraction and desorption temperature), a total of 26 experimental runs were performed, including six assays in the center of the cubic domain (pattern 00000), resulting in a total of 32 experiments. The CCD factor values and experiments are summarized in Supplementary Tables S1 and S2. The experimental data were fitted to a second-order polynomial equation and the model coefficients were determined based on a least-square regression analysis. The equations and the model suitability were assessed using the ANOVA approach and the main results are presented in Supplementary Table S3 (only the factors that significantly affect the response are shown). The R 2 values ranged from 0.80 to 0.97, which indicates a reasonable relationship between the experimental data and the fitted model. Models for most compounds show a F-probability below 0.1, indicating that variations in the response are associated with the models, rather than with the experimental error. However, M4Q and L5 show a higher F-probability, indicating a strong influence of experimental errors, possibly due to some problems related to competition within the fiber pores. Supplementary Table S4 presents the relevant variables and interactions identified by the Student's t-test, where the values in bold indicate the variables that present significant effect. Figure 2 presents the results from the Student's t-test for the main and quadratic effects. The main factors (X1 to X5) are the same defined in the screening design. None of the interactions presented Prob > |t| lower than 0.05. As can be observed, the "extraction time" (X2) and "extraction temperature" (X4) are the variables that significantly affect a larger number of responses, as defined by the low Student's t-test values.
of experimental errors, possibly due to some problems related to competition within the fiber pores. Supplementary Table S4 presents the relevant variables and interactions identified by the Student's t-test, where the values in bold indicate the variables that present significant effect. Figure 2 presents the results from the Student's t-test for the main and quadratic effects. The main factors (X1 to X5) are the same defined in the screening design. None of the interactions presented Prob > |t| lower than 0.05. As can be observed, the "extraction time" (X2) and "extraction temperature" (X4) are the variables that significantly affect a larger number of responses, as defined by the low Student's t-test values. Due to the high number of responses, the desirability function was applied to obtain the best compromise value for all variables, maximizing the response areas (Supplementary Figure S1). The optimal conditions were: 19.5% NaCl, 39 min extraction time, 10 min desorption time, 33 °C extraction temperature and 240 °C desorption temperature (D = 0.71). The optimization indicates that the addition of NaCl favors the transfer of analytes to the headspace and, consequently, to the fiber, due to the decrease of VMSs solubility in water (although it may also favor the co-extraction of interferents). Extraction time affects the amount of analyte extracted in the fiber coating, and this extraction is maximal when equilibrium is reached, achieving the highest sensitivity. However, in equilibrium extraction, competition might occur at higher concentrations. In this study, the optimal extraction temperature was set at 33 °C. Although the adsorption of fewer volatile compounds is favored by the increase in temperature (due to the improved mass transfer process of VMSs from sample to headspace), the partition coefficient between the coating and the sample decreases with increasing temperature, negatively affecting the adsorption of more volatile analytes. The desorption temperature of 240 °C is in the temperature range supported by the fiber and is high enough to properly desorb previously extracted molecules without damaging the fiber. Figure 3 presents an example of three-dimensional response surface plot for D6. Figure 3A presents the relationship between ionic strength and extraction time, Figure 3B the relationship between extraction time and desorption time, and Figure 3C the relationship between ionic strength and extraction temperature. As can be seen, the combined effect of the addition of salt (19.5% NaCl), together with the increase of desorption time (10 min) and extraction time (to 39 min), maximize the response for this compound. Due to the high number of responses, the desirability function was applied to obtain the best compromise value for all variables, maximizing the response areas (Supplementary Figure S1). The optimal conditions were: 19.5% NaCl, 39 min extraction time, 10 min desorption time, 33 • C extraction temperature and 240 • C desorption temperature (D = 0.71). The optimization indicates that the addition of NaCl favors the transfer of analytes to the headspace and, consequently, to the fiber, due to the decrease of VMSs solubility in water (although it may also favor the co-extraction of interferents). Extraction time affects the amount of analyte extracted in the fiber coating, and this extraction is maximal when equilibrium is reached, achieving the highest sensitivity. However, in equilibrium extraction, competition might occur at higher concentrations. In this study, the optimal extraction temperature was set at 33 • C. Although the adsorption of fewer volatile compounds is favored by the increase in temperature (due to the improved mass transfer process of VMSs from sample to headspace), the partition coefficient between the coating and the sample decreases with increasing temperature, negatively affecting the adsorption of more volatile analytes. The desorption temperature of 240 • C is in the temperature range supported by the fiber and is high enough to properly desorb previously extracted molecules without damaging the fiber. Figure 3 presents an example of three-dimensional response surface plot for D6. Figure 3A presents the relationship between ionic strength and extraction time, Figure 3B the relationship between extraction time and desorption time, and Figure 3C the relationship between ionic strength and extraction temperature. As can be seen, the combined effect of the addition of salt (19.5% NaCl), together with the increase of desorption time (10 min) and extraction time (to 39 min), maximize the response for this compound. Molecules 2021, 26, x 6 of 13 To prove the applicability of the empirical model, five additional tests with random patterns of the coded levels were performed. Results obtained were in the range of the predicted values by the proposed model and are presented in Supplementary Table S5.

Method Validation
Using the optimal conditions, the HS-SPME method performance was evaluated in terms of linearity, method detectability, precision, and accuracy. Tables 3-5 present the calculated validation parameters. Table 3. Linearity range and limits of detection and quantification for VMSs analysis by HS-SPME/GC-FID.

Analyte
Linearity Range (µg/L) (n = 8)  To prove the applicability of the empirical model, five additional tests with random patterns of the coded levels were performed. Results obtained were in the range of the predicted values by the proposed model and are presented in Supplementary Table S5.

Method Validation
Using the optimal conditions, the HS-SPME method performance was evaluated in terms of linearity, method detectability, precision, and accuracy. Tables 3-5 present the calculated validation parameters. Table 3. Linearity range and limits of detection and quantification for VMSs analysis by HS-SPME/GC-FID.

Analyte
Linearity Range (µg/L) (n = 8)  Regarding linearity of the response within the range of concentrations presented in Table 3, the correlation factor (R) was higher than 0.995 for all VMSs, except D3 and D6 and the relative standard deviation of the slope was less than 5%. Two different methods were used to evaluate the limits of detection. For linear siloxanes, the limits of detection (LODs) and quantification (LOQs) were estimated based on a signal-to-noise ratio (S/N) of 3 and 10, respectively. For cyclic siloxanes, they were determined based on 10 lab blanks, i.e., as the concentration of analyte that provides a response equal to the three and ten times their standard deviation, respectively. Nevertheless, these detection limits are still adequate for the intended screening approach using GC-FID.
Intra-day precision (repeatability) was evaluated through relative standard deviation (%RSD) of four replicates at different spiking levels (1 µg/L and 5 µg/L), while inter-day precision (intermediate precision) was assessed by determinations of the same sample, analyzed in different days (three days). In both cases, most relative standard deviation values were below 15%, as presented in Table 4. The mean relative standard deviation (RSD) values were about 12% for intra-day precision and about 15% for the inter-day precision, which is acceptable for this type for HS-SPME analysis.
Wastewater, tap water, and river water were used for accuracy measurements, which expresses the proximity of the obtained analytical result to the expected value. A standardaddition approach was used, by adding 25 ng of siloxanes to the samples and evaluating the resulting percentage of relatively recovery (%R defined as the ratio between the measured and the expected mass after the standard addition). Therefore, the concentrations of VMSs in the sample were calculated by the standard addition quantification method [14]. The matrix effect was previously reported for complex samples as wastewater, resulting in a signal suppression of 12-41% [12]. This was also verified in the present study. Therefore, to overcome this problem, a matrix-matched calibration was implemented for the analysis of wastewater samples. The obtained recoveries are also presented in Table 4. Examples of chromatograms may be seen in Supplementary Figure S2.
Finally, these parameters of validation were used to calculate the expanded uncertainty associated to final results accordingly to the protocol published by EURACHEM/CITAC Guide [15] and also Konieczka, and Namieśnik (2010) [16]. The global uncertainty was determined by identifying, estimating and combining all sources of uncertainty associated with the result [17,18]. The four main sources of uncertainty considered were: the uncertainty associated to the standards preparation (U1), calibration curve (U2), precision (U3) and accuracy (U4). Figure 4a shows the global uncertainty of D5, the compound referred in the literature as the most abundant in water. When concentrations decrease, approaching the limits of detection, the global uncertainty rises significantly. In Figure 4b the variation of the relative weight of each individual source may also be seen for D5. For the concentrations between 0.125 and 1 µg/L, U2 accounts for more than 50% of the global uncertainty. The standard preparation (U1) contributes around 10% or less for the global uncertainty. In the case of precision (U3) and accuracy (U4), they are increasingly higher as the upper concentration levels are reached, with their joint contribution reaching up to 70%. The global uncertainty for the other siloxanes follows the same trend as the presented in both figures (data not shown). To ensure a level of confidence of approximately 95%, results

Application of the Developed Method to Environmental Samples
Water samples were analyzed using the developed HS-SPME method. Results of the analyzed waters and respective uncertainty are presented in Table 5.
Siloxanes were not detected in both tap water and river water. In the wastewater sample, cyclic siloxanes presented higher concentrations, reaching up to 0.70 µg/L, which is in accordance with literature values. Linear siloxanes were detected with concentrations up to 0.44 µg/L. L3 and L4 concentrations were higher than those found in literature, while L5 is in accordance with the reported values, reaching up to 0.27 µg/L in influent [12]. The wastewater sample should be further analyzed by GC-MS, following the exact same procedure, in order to confirm the results.
In short, this methodology is suitable for the detection of VMSs in higher concentrations, such as effluents of industries that produce siloxanes or that utilize them in their working process. In addition, this methodology is particularly interesting for aqueous samples with low amounts of dissolved/suspended organic matter.
Individual stock solutions of each siloxane, including the internal standard M4Q, were prepared in n-hexane at 1 g/L. From those individual stock solutions, a 5 mg/L mix stock solution containing all the target analytes were prepared in acetone. A diluted M4Q individual stock solution, with a final concentration of 5 mg/L, was also prepared in acetone. From those stock solutions, new mix solutions with 50 µg/L and 500 µg/L were prepared in acetone. Eight calibration standards in acetone, with concentrations of each analyte ranging from 0.125 µg/L to 5 µg/L (internal standard concentration of 3 µg/L) were also prepared. All the solutions were stored protected from light and at −22 °C.

Application of the Developed Method to Environmental Samples
Water samples were analyzed using the developed HS-SPME method. Results of the analyzed waters and respective uncertainty are presented in Table 5.
Siloxanes were not detected in both tap water and river water. In the wastewater sample, cyclic siloxanes presented higher concentrations, reaching up to 0.70 µg/L, which is in accordance with literature values. Linear siloxanes were detected with concentrations up to 0.44 µg/L. L3 and L4 concentrations were higher than those found in literature, while L5 is in accordance with the reported values, reaching up to 0.27 µg/L in influent [12]. The wastewater sample should be further analyzed by GC-MS, following the exact same procedure, in order to confirm the results.
In short, this methodology is suitable for the detection of VMSs in higher concentrations, such as effluents of industries that produce siloxanes or that utilize them in their working process. In addition, this methodology is particularly interesting for aqueous samples with low amounts of dissolved/suspended organic matter.
Individual stock solutions of each siloxane, including the internal standard M4Q, were prepared in n-hexane at 1 g/L. From those individual stock solutions, a 5 mg/L mix stock solution containing all the target analytes were prepared in acetone. A diluted M4Q individual stock solution, with a final concentration of 5 mg/L, was also prepared in acetone. From those stock solutions, new mix solutions with 50 µg/L and 500 µg/L were prepared in acetone. Eight calibration standards in acetone, with concentrations of each analyte ranging from 0.125 µg/L to 5 µg/L (internal standard concentration of 3 µg/L) were also prepared. All the solutions were stored protected from light and at −22 • C.

Solid-Phase Microextraction Procedure
The starting point to define the operational variables was the published works that analyzed L3-L4 and D3-D6 in wastewater, and L2-L5 and D3-D6 in river water, respectively [12,13].
Ten mL water sample (containing 19.5% NaCl) was placed in a 20 mL screw cap amber vial (75.5 × 22.5 mm) with 18 mm, steel magnetic cap containing polytetrafluoroethylene (PTFE)/butyl septa. A small PTFE-coated stir bar was placed inside the vial and 30 ng of the internal standard (M4Q) was then added. A Supelco SPME holder for manual injection was used. Before the headspace solid-phase microextraction (HS-SPME) analysis, the sample vial was conditioned for 10 min in a bathwater, that was placed on top of the heated magnetic stirrer, at the extraction temperature of 33 ± 5 • C. Then, samples and calibration solutions were extracted with the SPME fiber at 33 • C for 39 min using a constant magnetic agitation rate of 750 rpm. Headspace mode of extraction was selected because VMSs are semi-volatile and due to the possible damage that the fiber could take during wastewater analysis, owing to the sample agitation and particulate matter. Finally, thermal desorption of the analytes was carried out by exposing the fiber in the GC injector port at 240 • C for 10 min.

GC-FID Analysis
The determination of the VMSs was carried out on a GC 2010 Plus gas chromatograph (Shimadzu, Kyoto, Japan). The chromatographic separation of the target compounds was performed on a Low-bleed DB-5MS ultra-inert column (5% phenyl, 95% methyl polysiloxane) fused silica capillary column (J&W Scientific, Folson, CA, USA), 30 m × 0.25 mm I.D., 0.25 µm film thickness. The oven temperature was programmed as follows: 35 • C hold for 5 min, raised at 10 • C/min to 95 • C, then 5 • C/min to 140 • C, and then 35 • C/min to 300 • C (hold for 5.5 min). The total time of analysis was 30 min. Helium was used as a carrier gas at a constant flowrate of 1 mL/min held by electronic flow control. The injector temperature was maintained at 240 • C and the splitless injection mode was used. The detector temperature was set at 300 • C. An SPME glass inlet liner (I.D., 0.75 mm, Shimadzu, Kyoto, Japan) and premium low bleed plasma coated septa (Shimadzu, Kyoto, Japan) were used.

Quality Assurance and Control (QA/QC)
Due to the VMSs ubiquity, analysts in the laboratory did not use personal care products and switched gloves whenever handling different samples. Glassware was subject to a special cleaning and decontamination procedure by pre-rinsing with acetone and distilled water. Non-calibrated material proceeded to further heating at 400 • C, for at least 1 h, in a laboratory muffle furnace. Procedural blanks (extraction of distilled water spiked with internal standard) were analyzed with every extraction batch, to correct sample concentration and identify any background contamination.

Design of Experiments
A design of experiments (DoE) consists of a series of runs that help to investigate the effects of input variables (factors) on an output variable (response). The factors that can influence the response are selected and evaluated. The main goal can be to identify significant factors that affect the response, identify factor interactions, to find optimal factor settings, match a target value or to build a predictive model, using the minimum number of experiments possible [19]. DoE is usually divided in three stages, starting with a screening design (identify the most important factors and reduce the number of factors that will be optimized), followed by a response surface design (construct a design that models a quadratic function of continuous factors, identifying factor settings that optimize the response) and finally the method validation.
Several models can be used for screening (e.g., full factorial and fractional factorial). Factorial designs are used to study the effects of two or more variables on the response.
Since the size of full factorial designs increases exponentially with the number of input variables, fractional factorial designs are preferred when dealing with high-dimensional problems. The performed 2 7-3 fractional factorial screening design (16 experiments, 1/8 of the full factorial design with 128 experiments) corresponds to a resolution IV design, meaning that the main effects are not confounded with other two-factor interactions. Another option would be to apply a 2 7−4 fractional factorial, that would result in only 8 experiments. However, this design corresponds to a resolution III design, which means it is possible to estimate the main effects, but these may be confounded with two-factor interactions. So, the 2 7-3 fractional factorial was preferred, even though it meant to conduct double the experiments. After performing the screening design and identifying the significant factors, a response surface methodology (RSM) is applied in order to model a quadratic function of continuous factors, identifying factor settings that optimize the response. Some examples include central composite design (CCD), Plackett-Burman design (PBD), Box-Behnken design (BBD), and Doehlert matrix design. The performed central composite designs consist of two-level full or fractional factorial design, coupled with two-star points (−α, +α) that represent new extreme values (low and high) for each factor in the design that allow the estimation of curvature. It also includes center points in order to measure stability and inherent variability of the model.
Each factor value (X i ) must be transformed into coordinates within a scale with dimensionless values (x i ), which must be proportional to its location in the experimental space. It enables the investigation of variables of different orders of magnitude without influencing the evaluation of the smallest. This means that the models operate on coded dimensionless input values x i (−1, +1) instead of actual factor values (X i ), that are converted by the following equation [20] (Equation (1)): where X 0 refers to the value of variable i in the center of the domain (x i = 0) and ∆X refers to the difference of that variable between x i = +1 and x i = 0. After conducting the experiments, the experimental responses must be adjusted to a quadratic model [20,21] (Equation (2)): where Y corresponds to the process response, x i to the codified independent variable, b 0 is the interception term, b i is the influence of the variable i in the response, b ii is the parameter that determines the shape of the curve, and b ij corresponds to the effect of the interaction among variables i and j. The quality of the fitted model can be evaluated by applying the analysis of variance (ANOVA), which is used to compare the variations due to the treatment with the variations due to random errors inherent to the measurements of the generated responses. Variations that occur in the response can be attributed to the model and are not due to random errors if the F-probability is less than 0.05 (for 95% confidence level). The relevant variables and interactions can be identified by the Student's t-test. If t-probability is smaller than 0.05, the parameter or interaction is considered significant. When several responses need to be evaluated (in this case, the main goal was to maximize the response areas for all compounds), a desirability function for each individual response is created, transforming each response into a dimensionless individual desirability (d i ) scale, that ranges from d i = 0 for inacceptable value to d i = 1 for the single response maximum. The response is considered unacceptable if it is outside the desirable limit [20] (Equation (3)): where d i is the individual desirability values, f (x) the response value, A is the minimum response value, B the maximum response value, and w is the weight used to determine the importance (in this work, equal weights were chosen for all responses). It is then possible to calculate the overall desirability (D), also ranging from 0 to 1 [21] (Equation (4)): where m is number of responses studied in the optimization process. The optimized conditions indicate the values of input variables for which the overall desirability is maximal. In order to generate the experimental matrix for the design of experiments, to evaluate the effects of several parameters and to perform the data analysis, the JMP 14.3 (SAS Institute Inc., Cary, NC, USA) statistical software was used. However, this could be substituted by other similar statistical software.

Conclusions
This work demonstrates the relevance of experimental design and method validation. Indeed, it goes further, by linking the two components and applying them to the analysis of VMSs in water. The application of design of experiments (DOE) enabled the study of the effect of input variables (factors) on an output variable (response) using the minimum number of experiments possible, thus saving time and resources. The developed HS-SPME technique emerges as an interesting solution to conventional procedures since it combines the potential to eliminate solvent consumption and the concomitant issues of used solvent disposal, with a simple, low-cost, and eco-friendly extraction procedure. In terms of performance parameters, the developed HS-SPME method presents similar results to those obtained by other studies found in the literature that also used HS-SPME for the analysis of VMSs in water samples. The present work also has the advantage of using a simple and cheaper instrumental methodology based on GC-FID, which is more easily replicable in WWTP analysis laboratories, to screen a larger number of samples. Table S1: Factors and values for the proposed central composite design); Table S2: Proposed central composite design); Table S3: Second-order polynomial equation obtained for each target compound and model suitability parameters for the response functions; Table S4: Results from the Student's t-test for the main, quadratic and interaction effects; Table S5: Additional testing to prove the applicability of the proposed empirical model; Figure S1: Optimal conditions and desirability function; Figure

Data Availability Statement:
The data presented in this study is available on request from the corresponding author.