Computational-Based Study of QuEChERS Extraction of Cyclohexanedione Herbicide Residues in Soil by Chemometric Modeling

Assessment of two buffered QuEChERS (quick, easy, cheap, effective, rugged, and safe) versions (i.e., citrate and acetate) modified by including methanol to recover the residues of three cyclohexanedione oxime (CHD) herbicides and three of their byproducts from agricultural soil was performed. In this context, a full second-order face-centered factorial experimental design was developed to quantify the influences of the main five variables (i.e., extraction time, water content, soil weight, and extraction solvent volume and composition) on the target compound recoveries. The fitting equations satisfactorily described the extraction process behavior. The mathematical models also showed the most influencing independent variables (i.e., extraction solvent composition and soil weight). Handling simpler expressions was possible with the acetate QuEChERS but not with the citrate QuEChERS. The recoveries of the CHD residues were close to 100% after performing the extraction under suitable conditions. Furthermore, dispersive solid-phase extraction (dSPE) clean-up steps were assessed to reduce the matrix effect in mass spectrometry. In this sense, the citrate QuEChERS in combination with the PSA + C18 clean-up step was the best option for the extraction of CHD residues.


Introduction
Pesticides are ubiquitous contaminants in the environment [1][2][3], as field soils are a frequent receptor of these compounds when they are applied to crops [4,5]. Not only pesticides but also their degradation products-which result from the environmental exposure of pesticides to microorganisms, light, etc.-are found in soil. Therefore, the establishment of reliable methods for analyzing pesticide residues in soils is of utmost importance to assess their effects on this environmental compartment [3,6,7]. However, the determination of pesticide residues in soils is not a simple task due to the complexity and heterogeneity of this matrix as well as to the low concentrations of pesticide residues commonly found.
Cyclohexanedione oxime (CHD) herbicides are thermolabile and low volatility pesticides used for postemergence control of annual and perennial grass weeds in broad-leaved crops such as sugar beet, soybean, or oilseed rape. In general, CHD herbicides are rapidly degraded under different environmental conditions, including biotic and abiotic processes, leading to the formation of new byproducts [8][9][10]. These new compounds may pose a potential threat to succeeding crops or nontarget organisms. Indeed, different studies have suggested the possibility that some of the byproducts Figure 1. Chemical structure of the selected CHD herbicides (i.e., alloxydim, sethoxydim and profoxydim) and three of their degradation products (i.e., deallyloxylated-alloxydim, deethoxylatedsethoxydim and sethoxydim-oxazole).

Optimization of the Chromatographic Method
First, the chromatographic conditions were optimized to achieve the optimum separation of the six target compounds.
Due to the acidic character of the target compounds, the addition of formic acid helps improve the protonation of the target compounds and hence the peak shape. Thus, the best results were obtained using a mixture of acetonitrile-water in which 0.1% formic acid was used to acidify the mobile phase. The addition of a higher concentration of formic acid did not improve the separation or mass ionization.
Regarding the mass analyzer, the best mass response for all target compounds was observed in Electrospray Ionization (ESI) positive ion mode. Moreover, the most intense ions for the target compounds corresponded to their protonated molecule ions [M + H] + , so these ions were selected for quantification purposes/matrix study (Table S4).

Modeling of the QuEChERS Extraction Method
For the development of the factorial design, photodiode array detector (DAD) was selected as the detector due to the variability of the matrix effect on the mass response observed under the different extraction conditions. Thus, the use of mass spectrometry (MS) detectors would make it difficult to model the QuEChERS method since it would involve the preparation of 90 different external matrix-matched calibrations.
The aim of the factorial design is to identify the key variables and possible synergistic effects controlling the process [21,23], which allow modeling of the QuEChERS extraction to maximize the recoveries of CHD residues from soil. Bruzzoniti et al. [24] have already highlighted the necessity of using chemometric tools to contribute to a better understanding of the extraction process in soil.
Tables S1 and S2 show the recoveries obtained for the CHD residues during the factorial design of the QuEChERS extraction. The recoveries for the majority of compounds are usually within the  Figure 1. Chemical structure of the selected CHD herbicides (i.e., alloxydim, sethoxydim and profoxydim) and three of their degradation products (i.e., deallyloxylated-alloxydim, deethoxylatedsethoxydim and sethoxydim-oxazole).

Optimization of the Chromatographic Method
First, the chromatographic conditions were optimized to achieve the optimum separation of the six target compounds.
Due to the acidic character of the target compounds, the addition of formic acid helps improve the protonation of the target compounds and hence the peak shape. Thus, the best results were obtained using a mixture of acetonitrile-water in which 0.1% formic acid was used to acidify the mobile phase. The addition of a higher concentration of formic acid did not improve the separation or mass ionization.
Regarding the mass analyzer, the best mass response for all target compounds was observed in Electrospray Ionization (ESI) positive ion mode. Moreover, the most intense ions for the target compounds corresponded to their protonated molecule ions [M + H] + , so these ions were selected for quantification purposes/matrix study (Table S4).

Modeling of the QuEChERS Extraction Method
For the development of the factorial design, photodiode array detector (DAD) was selected as the detector due to the variability of the matrix effect on the mass response observed under the different extraction conditions. Thus, the use of mass spectrometry (MS) detectors would make it difficult to model the QuEChERS method since it would involve the preparation of 90 different external matrix-matched calibrations.
The aim of the factorial design is to identify the key variables and possible synergistic effects controlling the process [21,23], which allow modeling of the QuEChERS extraction to maximize the recoveries of CHD residues from soil. Bruzzoniti et al. [24] have already highlighted the necessity of using chemometric tools to contribute to a better understanding of the extraction process in soil.
Tables S1 and S2 show the recoveries obtained for the CHD residues during the factorial design of the QuEChERS extraction. The recoveries for the majority of compounds are usually within the acceptable range of 70-110% at the EU level for analytical methods according to Regulation (EU) No 283/2013 (Figure 2, Tables S1 and S2). An exception to this generality is observed in several citrate tests for sethoxydim and profoxydim for which the recoveries surpass the regulated threshold. These observations seem to be in agreement with those made by Lehotay et al. [25], who found that the acetate-buffered version of QuEChERS gives better recoveries for most of the pesticides studied.  Tables S1 and S2). An exception to this generality is observed in several citrate tests for sethoxydim and profoxydim for which the recoveries surpass the regulated threshold. These observations seem to be in agreement with those made by Lehotay et al. [25], who found that the acetate-buffered version of QuEChERS gives better recoveries for most of the pesticides studied.

Figure 2.
Experimental results and correlation between citrate and acetate quick, easy, cheap, effective, rugged, and safe (QuEChERS) recoveries for each of the extraction conditions assessed. Figure 2 shows a large dispersion regarding the recoveries of the target compounds. Therefore, precise selection of the extraction conditions is required to obtain suitable recoveries. Equation 1 fitted the experimental data well by means of least squares multiple regressions, and the resultant polynomial models are shown in Tables 1 and 2 together with the significance levels of the factors and the statistics parameter displaying the goodness of fit.

Factorial Design of Citrate QuEChERS
In the case of citrate QuEChERS, AM is the variable that exerts the largest influence on recoveries at a confidence level of 95% (Table 1). For all target compounds except for sethoxydim-oxazole, when AM increases, the recoveries decrease (Table 1), especially in the high AM range ( Figure 3). However, a high value of WC can change this trend for CHD herbicides, as shown in Figure 4a, where the typical shape of a response surface for the target compounds is depicted. The negative sign of the fitting parameter bWC for these CHD herbicides (Table 1) is responsible for the response surface falling in the area of high WC. This fact seems opposite to the assumed importance of the water to the successful extraction of the pesticides from the soil [19,[26][27][28]. However, bWC is significant only in the best case at a confidence level of 80% (Table 1), and these observations may be related to the inherent error associated with the experimental procedure. In this sense, the decrease in recoveries can be easily offset by increasing the AM slightly ( Figure 4a). In the case of sethoxydim-oxazole, a significant opposite trend is observed for AM in whole range of WC with a comparatively higher intensity ( Table  1, Figures 3 and 4). The comparison performed in Figure 4 between sethoxydim-oxazole and, for example, profoxydim, clearly shows the highlighted differences between sethoxydim-oxazole and the rest of the compounds.   Tables 1 and 2 together with the significance levels of the factors and the statistics parameter displaying the goodness of fit.

Factorial Design of Citrate QuEChERS
In the case of citrate QuEChERS, AM is the variable that exerts the largest influence on recoveries at a confidence level of 95% (Table 1). For all target compounds except for sethoxydim-oxazole, when AM increases, the recoveries decrease (Table 1), especially in the high AM range ( Figure 3). However, a high value of WC can change this trend for CHD herbicides, as shown in Figure 4a, where the typical shape of a response surface for the target compounds is depicted. The negative sign of the fitting parameter b WC for these CHD herbicides (Table 1) is responsible for the response surface falling in the area of high WC. This fact seems opposite to the assumed importance of the water to the successful extraction of the pesticides from the soil [19,[26][27][28]. However, b WC is significant only in the best case at a confidence level of 80% (Table 1), and these observations may be related to the inherent error associated with the experimental procedure. In this sense, the decrease in recoveries can be easily offset by increasing the AM slightly ( Figure 4a). In the case of sethoxydim-oxazole, a significant opposite trend is observed for AM in whole range of WC with a comparatively higher intensity ( Table 1, Figures 3 and 4). The comparison performed in Figure 4 between sethoxydim-oxazole and, for example, profoxydim, clearly shows the highlighted differences between sethoxydim-oxazole and the rest of the compounds.
Moreover, an important significant value at a 95% confidence level for sethoxydim-oxazole is also observed in the fitting parameter associated with the quadratic term of AM (b AM-AM ) with a negative value (Table 1), which causes the existence of maxima in the graphical representation of the AM and recovery values (Figures 3 and 4b). Therefore, good recoveries of sethoxydim-oxazole are achieved at high AM values, and its recoveries are very sensitive to small variations of the AM variable.                  On the other hand, two important interactions between AM and WC and between AM and EV (bAM-WC and bAM-EV, respectively) are observed on the target compound recoveries (Table 1). Thus, except for sethoxydim-oxazole, an important positive effect is observed for the fitting parameter bAM-WC (Table 1), which could be considered a synergistic effect on the recoveries of all target compounds due to the negative effect and high importance of the individual fitting parameter bAM and the low importance of the individual fitting parameter bWC in the polynomial model (Table 1). This effect of the bAM-WC fitting parameter on favoring the dissolution of the polar CHD residues might come from an increase in polarity and protic characteristics upon reaching a sufficient threshold in the extraction medium when acetonitrile, in combination with methanol, is added and combined with water [29,30]. Conversely, a negative effect on the target compound recoveries is observed for the fitting parameter bAM-EV (Table 1). However, in this case, the AM effect on the target compound recoveries is partially compensated for by the EV (i.e., acetonitrile and methanol) by adding more quantity of methanol, that is, a solvent that is polar and protic but less so than water [29,30], thereby diminishing the negative significance of the fitting parameter bAM-EV regarding the significance of the individual fitting parameter bAM. The fitting parameter bAM-EV seems to indicate a synergistic effect in the case of On the other hand, two important interactions between AM and WC and between AM and EV (b AM-WC and b AM-EV , respectively) are observed on the target compound recoveries (Table 1). Thus, except for sethoxydim-oxazole, an important positive effect is observed for the fitting parameter b AM-WC (Table 1), which could be considered a synergistic effect on the recoveries of all target compounds due to the negative effect and high importance of the individual fitting parameter b AM and the low importance of the individual fitting parameter b WC in the polynomial model (Table 1). This effect of the b AM-WC fitting parameter on favoring the dissolution of the polar CHD residues might come from an increase in polarity and protic characteristics upon reaching a sufficient threshold in the extraction medium when acetonitrile, in combination with methanol, is added and combined with water [29,30]. Conversely, a negative effect on the target compound recoveries is observed for the fitting parameter b AM-EV (Table 1). However, in this case, the AM effect on the target compound recoveries is partially compensated for by the EV (i.e., acetonitrile and methanol) by adding more quantity of methanol, that is, a solvent that is polar and protic but less so than water [29,30], thereby diminishing the negative significance of the fitting parameter b AM-EV regarding the significance of the individual fitting parameter b AM . The fitting parameter b AM-EV seems to indicate a synergistic effect in the case of sethoxydim-oxazole: both independent variables, i.e., AM and EV are significant at a 0.05 level with a positive effect on sethoxydim-oxazole recovery, while their interaction is also significant at a 0.05 level but with a negative effect in this case (Table 1). This effect might be observed when the medium becomes sufficiently polar and protic. The interaction of the three variables, i.e., AM, EV and WC (b AM-EV-WC ) ( Table 1), with negative effects on the recoveries of the CHD residues should mainly result from acetonitrile.

Factorial Design of Acetate QuEChERS
In this case, a significant global influence is given by SW at a confident level of 95% on the recoveries for all target compounds ( Table 2). The greater influence and tendency observed with the variable AM in the citrate QuEChERS appears to be maintained in the acetate version for all compounds except for the dealkoxylated byproducts (i.e., from alloxydim and sethoxydim). Regarding the WC, this variable has more influence on the degradation products of sethoxydim (i.e., deethoxylated-and oxazole-sethoxydim) at a confidence level of 95%, followed by profoxydim and deallyloxylated-alloxydim at a confidence level of 90%. Moreover, the WC shows a higher positive influence on recoveries of all target compounds in the acetate QuEChERS versus the citrate QuEChERS.   Figure 3). However, in this case, significance is observed in the AM quadratic term (b AM-AM ) for all target compounds ( Table 2). The EV quadratic terms (b EV-EV ) for alloxydim and sethoxydim are also statistically significant (0.05 and 0.10 level, respectively) but in this case with a positive sign (Table 2). Therefore, both AM and EV should be taken into special consideration during the extraction of sethoxydim and alloxydim because maxima and minima are observed when both variables are depicted together, resulting in maximal recoveries at both extremes of the EV tested and with medium or low AM (see an example in Figure 5). This effect is not so evident for the rest of the target compounds. In contrast to acetate, citrate QuEChERS did not exhibit this trend since all quadratic terms were negative (Table 1); i.e., there were no minima. In this sense, Figure 5 shows that once, for example, the maximum extraction at the lowest EV is reached, if the AM is kept constant while the EV is increased, then the greater polar and protic effect from the water of the medium relative to the lower polar and protic effect that a higher EV implies may "vanish" since methanol and acetonitrile are less polar and protic than water. However, the amount of methanol in the medium is eventually sufficient to compensate for the lower polar and protic power of the extraction medium; thus, the maximum extractions are reached again. This effect, the importance of the individual fitting parameter b WC (Table 2) and the synergistic effect favoring the target compound recoveries at a 0.05 level between SW and WC, which is observed when the individual fitting parameters b SW and b WC are compared with the fitting parameter b SW-WC (Table 2), clearly denote the importance of water in the case of the acetate QuEChERS. A plausible explanation for the synergistic effect among SW and WC might come from soil rehydration, which improves the access to the soil pores and, consequently, the pesticide residue partitioning process toward acetonitrile/methanol. Competition between water and some pesticides for adsorption sites of soil humic substances by forming hydrogen bonds can also be a simultaneous process that promotes pesticide desorption. These two hypotheses were also provided by other authors [27,31] after observing an improvement in the recoveries of pesticides with different QuEChERS procedures after soil hydration. the importance of the individual fitting parameter bWC (Table 2) and the synergistic effect favoring the target compound recoveries at a 0.05 level between SW and WC, which is observed when the individual fitting parameters bSW and bWC are compared with the fitting parameter bSW-WC (Table 2), clearly denote the importance of water in the case of the acetate QuEChERS. A plausible explanation for the synergistic effect among SW and WC might come from soil rehydration, which improves the access to the soil pores and, consequently, the pesticide residue partitioning process toward acetonitrile/methanol. Competition between water and some pesticides for adsorption sites of soil humic substances by forming hydrogen bonds can also be a simultaneous process that promotes pesticide desorption. These two hypotheses were also provided by other authors [27,31] after observing an improvement in the recoveries of pesticides with different QuEChERS procedures after soil hydration. A high ET when the AM is increased has shown detrimental effects on the recoveries of the target compounds (see fitting parameter bAM-ET in Table 2), possibly as a consequence of increased diffusion and adsorption of the CHD residues in/on solids, i.e., soil and salts of acetate QuEChERS.
Finally, the different behavior of sethoxydim-oxazole in both extraction QuEChERS (i.e., citrate and acetate) seems to indicate that the absence of a keto-enol equilibrium in sethoxydim-oxazole, which conversely occurs in the rest of target compounds (see chemical structures in Table S4), influences its interaction either with the soil or with the solvents. A high ET when the AM is increased has shown detrimental effects on the recoveries of the target compounds (see fitting parameter b AM-ET in Table 2), possibly as a consequence of increased diffusion and adsorption of the CHD residues in/on solids, i.e., soil and salts of acetate QuEChERS.
Finally, the different behavior of sethoxydim-oxazole in both extraction QuEChERS (i.e., citrate and acetate) seems to indicate that the absence of a keto-enol equilibrium in sethoxydim-oxazole, which conversely occurs in the rest of target compounds (see chemical structures in Table S4), influences its interaction either with the soil or with the solvents.

Handling of the Polynomial Models
To handle simpler expressions, the nonsignificant terms (e.g., below the 95% confidence level) of the obtained empirical mathematical equations could be removed [32]. However, a proper demonstration of the suitability of this action is required. In this sense, to avoid unreal distortions as a consequence of an improper definition of the b 0 constant, a study to assess the effect of removing terms according to their significance level from the empirical correlations was planned. Table 3 shows the existing deviations between the experimental recoveries and the polynomial model results determined according to the minimum fixed confidence level established (i.e., confidence level ≥95%, ≥90%, ≥80%, ≥70% and ≥0% (full equation)).  The deviations obtained after removing all terms of the polynomial models below the 95% confidence level are always below 10% for the acetate QuEChERS. However, this action is not suitable in the case of citrate QuEChERS. In fact, at the 70% confidence level, none of the six target compounds are determined with acceptable deviations within any of the ranges studied for the different independent variables (Table 3). Moreover, the deviations between the experimental recoveries and the polynomial model results were lower for the acetate QuEChERS than for the citrate QuEChERS before and especially after removing the nonsignificant terms. The lower citrate response surface curvature, defined by less significant quadratic and interaction terms (Tables 1 and 2), explains the differences obtained with both QuEChERS. These results clearly indicate that a case-by-case study should be performed if the intention is to remove terms from the polynomial models. However, keeping all terms in the equations does not make the model incorrect; instead, it makes the model more complex but also more accurate.
On the other hand, if the R values showing the goodness of fit are compared, it is clear that the R values are not in agreement with the range of deviations. In this sense, a low R in the QuEChERS acetate does not mean a great deviation, e.g., deallyloxylated-alloxydim and deethoxylated-sethoxydim at the 95% confidence level. This is due to a low dispersion in the recovery values, with response surfaces that cross the whole cloud of recovery points in a suitable way, achieving low deviations relative to the experimental results. An opposite situation is observed, for example, with sethoxydim-oxazole in the citrate QuEChERS at any confidence level (Table 3). In this case, the wide range in the recovery values provides a high R, while the presence of some individual recoveries far away from the response surfaces is responsible for the worst deviations [33]. Therefore, the regression parameter R alone should not be considered sufficient to establish the accuracy of the models.

Confirmation of Extraction Conditions and Clean-Up Assessment
According to the experimental results previously discussed, the polynomial models obtained for citrate and acetate QuEChERS (Tables 1 and 2) provide different extraction conditions to achieve recoveries of 100% for the target compounds when full equations are used. Therefore, the best experimental management conditions were selected for the extraction. In this sense, as an additional criterion, the minimum values of the variables were selected for the confirmation studies. Table 4 shows the final optimized extraction conditions determined by the polynomial models and the experimental recoveries obtained. Despite advances in sample preparation methods, pesticide residues are often found at levels too low to be detected by DAD detectors in environmental matrices. In this regard, mass spectrometers are known to provide a higher sensitivity and selectivity, but the matrix effect is a major drawback for the quantitative determination of analytes by this technique [3,34]. Thus, once the best conditions for citrate and acetate QuEChERS extraction were selected (Table 4), the matrix effect was evaluated under these conditions to select the best combination of QuEChERS extraction/clean-up step (PSA and PSA-C18) to analyze CHD residues by HPLC-MS.
As illustrated in Figure 6, suppressive effects on MS signals are observed for all the target compounds when either the QuEChERS citrate or acetate extraction kit is used without a clean-up step. This signal suppression is higher with citrate QuEChERS than with the acetate version. Although the best results were observed with the acetate QuEChERS, the coelution of soil components makes it necessary to subtract the matrix blank signal for quantification and to thoroughly clean the chromatographic column after each injection. Despite advances in sample preparation methods, pesticide residues are often found at levels too low to be detected by DAD detectors in environmental matrices. In this regard, mass spectrometers are known to provide a higher sensitivity and selectivity, but the matrix effect is a major drawback for the quantitative determination of analytes by this technique [3,34]. Thus, once the best conditions for citrate and acetate QuEChERS extraction were selected (Table 4), the matrix effect was evaluated under these conditions to select the best combination of QuEChERS extraction/clean-up step (PSA and PSA-C18) to analyze CHD residues by HPLC-MS.
As illustrated in Figure 6, suppressive effects on MS signals are observed for all the target compounds when either the QuEChERS citrate or acetate extraction kit is used without a clean-up step. This signal suppression is higher with citrate QuEChERS than with the acetate version. Although the best results were observed with the acetate QuEChERS, the coelution of soil components makes it necessary to subtract the matrix blank signal for quantification and to thoroughly clean the chromatographic column after each injection.  To reduce the matrix effect, two different clean-up QuEChERS kits were used, PSA and PSA-C18. PSA removes sugars and fatty acids, whereas the C18 sorbent effectively traps and removes lipids, starch, and humic substances [19,35,36]. In this study, two opposite results were observed: whereas citrate QuEChERS extraction combined with PSA or PSA-C18 clean-up led to a reduction in the suppressive matrix effect on all the target compounds (Figure 6), the combination of acetate with both clean-up QuEChERS led to a higher suppression of the MS signals. Notably, with the citrate QuEChERS, a similar reduction in the matrix effect was observed for both clean-up PSA and PSA-C18. Regarding the acetate QuEChERS, the coelution of soil components was reduced when PSA-C18 was used as clean-up step.
According to the findings mentioned above, the best combination for the determination of CHD residues from soil entails extraction with citrate QuEChERS followed by the PSA-C18 clean-up step. However, since the matrix effect is maintained above 20%, the use of external matrix-matched standards is recommended [37,38].
Acetonitrile (HPLC superGRAD grade) and methanol (HPLC grade) were acquired from Macron Fine Chemicals (Gliwice, Polland). The water used for the LC mobile phase and the aqueous solutions were purified with a Millipore system (Milli-Q-50 18 mΩ) (Molsheim, France).
Stock solution containing a mixture of the six target compounds (i.e., alloxydim, deallyloxylatedalloxydim, sethoxydim, sethoxydim-oxazole, deethoxylated-sethoxydim and profoxydim) was prepared in methanol at a concentration of 500 mg L −1 . This solution remained stable for ten days when stored in the dark at −18 • C. Working solutions were prepared daily by dissolving the appropriate amount of the composite stock solution in acetonitrile/methanol or in matrix extracts.

Soil Samples
The soil used in the experiments was collected from an agricultural field located in Sevilla province (southwest of Spain) with no history of CHD herbicide application. The soil was taken from the top layer (0-20 cm), and after being air-dried at room temperature, it was passed through a 2 mm sieve and stored at 4 • C until use. Soil texture was determined by sedimentation using the pipette method [40]. The organic carbon an nitrogen content were determined according to the Walkley and Black method [41] and the Kjeldahl method [42], respectively. Soil pH was measured in a 1:2.5 (w/v) soil/deionized water mixture. Soil at the site was classified as sandy clay loam containing approximately 32% sand, 55% silt, and 13% clay, with organic carbon = 9.3 g kg −1 dw, N Kjeldahl = 0.89 g kg −1 dw and pH = 7.2.
The soil samples were spiked with the stock solution containing 10 mg kg −1 dw of each of the six target compounds, which corresponds to a field application rate of 200 mg ha −1 of each target compound, close to the recommended application doses for CHD herbicides. The mixture was vigorously shaken for 1 min with an Intelli-Mixer device (JP SELECTA S.A., Barcelona, Spain) and subsequently kept in the dark for 2 h to allow the interaction between the target compounds and the matrix.

QuEChERS Procedure
The general QuEChERS method was as follows: a portion of soil was weighed (SW) into a 50 mL centrifuge tube, and water was added (WC). Afterwards, different EVs of solvent mixtures (acetonitrile/methanol (AM)) were added, and the mixture was shaken for 1 min with the Intelli-Mixer device. The content of the commercial extraction QuEChERS (citrate or acetate) was added, and the mixture was immediately shaken for different ETs and subsequently centrifuged for 5 min at 4000 rpm and 4 • C.
For the clean-up step by dispersive solid-phase extraction (dSPE), an aliquot of the supernatant of the QuEChERS extraction tube (2 mL) was transferred into a 15 mL centrifuge tube containing commercial clean-up QuEChERS (PSA or PSA-C18). The tube was shaken for 1 min and centrifuged for 5 min at 4000 rpm and 4 • C, and the resulting supernatant was filtered using a 0.45 µm PVDF syringe filter for chromatographic injection.

Factorial Design of QuEChERS Extraction
The main independent variables selected during the QuEChERS extraction and their working ranges were the following: ET, (1-3 min.); WC, (1.5-4.5 mL); AM, (5-95%); EV, (8-16 mL) and SW, (2.4-5.8 g). The SW range was maintained below 6 g to facilitate the management of a representative soil sample in the lab, and the minimum EV was selected as a function of the wettability of the QuEChERS kit and the concentrations of the target compounds in the final extracts.
3 5 full second-order face-centered factorial experimental designs with three center points were performed to determine the influence of selected factors and their effects on the extraction of target compounds from soil [33].
The general model for recovery is the following polynomial equation: where Y represents the dependent variable or system response (recovery), b 0 is a constant that fixes the response at the central point of the experiment, bi is the fitting parameter for the linear effect term, and b ii , b ij and b ijk are the fitting parameters for the interaction effects terms. These fitting parameters allow the determination and comparison of the effects of each independent variable on recovery because of normalization. X i , X j and X k represent the normalized independent variables previously defined. Statistical calculations and analysis were performed using the Excel statistical module. Tables S1 and S2 summarize the structure of the experimental design case by case.

Chromatographic Analysis
Chromatographic analysis was performed with an HPLC system (series 1100; Agilent Technologies, Palo Alto, USA) coupled with a DAD and a single-quadrupole mass spectrometer equipped with an ESI interface. The separation was achieved at 20 • C on a Phenomenex Kinetex ® 2.6 µm C18 100 Å (100 × 4.6 mm) column protected by an AJO-4287 guard cartridge. Elution solvents were ultrapure water acidified with 0.1% formic acid (A) and acetonitrile (B), and the elution was performed at a flow rate of 0.7 mL min −1 using the following gradient method: 50% B at 0 min, increase 50-95% B in 3 min and hold at 95% for 4 min.
MS quantification was performed using the most intense product ion ( Figure S1), and the mass spectrometer operated in positive mode under the following conditions: nebulizer pressure: 35 psi; drying gas flow: 12.0 L min −1 ; drying gas temperature: 350 • C; gain: 1; nozzle voltage: 3000 V and fragmentor voltage: 150 V.
To evaluate the matrix effect in MS, the response of each of the six target compounds in the extraction solvent was compared with the response of these compounds in blank soil extracts within the linear concentration range (0.005-0.1 mg L −1 ). Blank soil extracts were obtained from the optimized extraction conditions with both QuEChERS extraction kits (i.e., citrate and acetate) alone and in combination with the clean-up steps tested (i.e., PSA and PSA-C18). The matrix effect on the target compounds was studied according to Equation 2 by comparing the slopes between the calibration curves in the extraction solvent and in the soil matrix. ME(%) = S m − S s S s · 100 (2) where ME is the matrix effect, S m is the slope in the soil matrix and S s is the slope in the extraction solvent. ME values below 0% indicate signal suppression, whereas values above 0% indicate signal enhancement. Notably, an ME > |20%| is commonly considered to have a significant impact on the performance of the analytical method.

Conclusions
A full 3 5 second-order face-centered design of experiments was performed to examine in detail the effects of five independent variables, i.e., ET, WC, EV, SW and AM, on two different QuEChERS extraction systems, i.e., citrate and acetate, to efficiently recover three CHD herbicides and three of their byproducts from agricultural soil. This study involved 90 different experimental conditions between both systems, and the empirical mathematical models were found to fit well with the experimental results, yielding the following main results: highly polar and protic extraction media favor the target compound recoveries, while the AM exerts a significant and negative influence. Sethoxydim-oxazole follows the opposite rules. In general, the influence of the independent variables follows the order AM > SW > WC ≥ EV >> ET. A proper set of independent variables was selected, which enables recoveries close to 100% for the six target compounds studied: SW = 3.23 g, AM = 52.8%, EV = 8 mL, WC = 3.8 mL and ET = 1.8 min for citrate QuEChERS and SW = 3.96 g, AM = 57.0%, EV = 9.4 mL, WC = 4.0 mL and ET = 2.8 min for acetate QuEChERS.
QuEChERS does not readily achieve good recoveries with highly polar pesticides with log K ow < −2 as a result of their poor or no partition into the organic phase. The aforementioned conclusions indicate that both QuEChERS versions, modified by introducing an additional solvent into the extraction system, i.e., methanol (polar/protic solvent), are good options for the simultaneous determination of pesticides with different behavior toward acetonitrile (polar/aprotic solvent) alone. In view of these promising results, an assessment of these QuEChERS versions over simultaneous residue recoveries of highly polar pesticides (e.g., fosetyl, ethephon and glyphosate) and nonpolar pesticides (e.g., pyrethroids and organophosphates such as chlorpyrifos) should be the next step toward the global use of these new QuEChERS approaches.
It is suggested that if nonsignificant terms want to be removed from the polynomial models, then a pertinent justification should be provided. In this particular case, the equation for acetate QuEChERS extraction can be simplified without great deviations relative to the experimental results. However, the results with the citrate QuEChERS suggest that all terms should be kept in the equation.
The extraction process was further studied by assessing the matrix effect observed when mass spectrometry instruments were used. Citrate QuEChERS in combination with the PSA + C18 clean-up step is more efficient in reducing the matrix effect than the acetate version with a clean-up step.
In summary, this methodology would allow optimizing the analytical procedure to analyze pesticide residues that possess different behavior during the extraction step. The application of chemometrics models can minimize analysis time, size of sample and solvent expense without compromising the analytic performance that will improve the accuracy of the method and thus the risk assessment of pesticides.
Supplementary Materials: Figure S1 and Tables S1-S4 were provided as supplementary material and are available online. Funding: Funding for this research was provided by the CICYT project RTA2014-00044-00-00 and RTA2017-00043-00-00.

Conflicts of Interest:
The authors declare no conflicts of interest.