Innovative Strategy for Aroma Stabilization Using Green Solvents: Supercritical CO2 Extracts of Satureja montana Dispersed in Deep Eutectic Solvents

The aim of this work was to establish the potential of natural deep eutectic solvents (NADES) for the stabilization of aroma volatile organic compounds from a natural source. Satureja montana was used as a source of volatile components, as it is rich in terpenes of great commercial and biological importance, such as carvacrol, thymol, and thymoquinone, among others. Supercritical CO2 was used to extract the lipophilic fraction of S. montana, which was further directly dispersed in NADES. The stabilizing capacity of seven different NADES based on betaine and glycerol was analyzed. The stability of the components in NADES was monitored by analyzing the headspace profile during 6 months of storage at room temperature. The changes in the headspace profile over time were analyzed by using different statistical and chemometric tools and the Wilcoxon matched pair test. It was determined that alterations over time occurred such as degradation and oxidation, and they were the most prominent in the control. In addition, the indicator of decreased stability of the control was the formation of the new compounds that could compromise the quality of the product. In the stabilized NADES samples, the changes were significantly less prominent, indicating that the NADES had a stabilizing effect on the volatile compounds. According to Wilcoxon matched pair test, the most efficient stability was achieved by using betaine/ethylene glycol, glycerol/glucose, and betaine/sorbitol/water. Therefore, by applying two green solvents, a sustainable approach for obtaining pure and high-quality S. montana extracts with extended stability at room temperature was established.


Introduction
Volatile organic aroma compounds represent important constituents in various fields because they possess sensory and numerous biologically important properties. These components are incorporated into products intended for the medical treatment of humans and animals, cosmetic and food products, beverages, fragrances and perfumes, products for care and cleaning, insecticides, repellents, etc. Several markets drive the demand of volatile organic aroma compounds, including the food and beverage (35%), fragrances, cosmetics and aromatherapy (29%), household (16%), and pharmaceutical (15%) markets [1,2]. Numerous industries use synthetic compounds, as the process of obtaining them from natural sources can be expensive and irrational. However, the connection of synthetic components with negative consequences on health and the environment, together with the growth of green consumerism, represent the reason for the growing global demand for natural aroma compounds, their implementation in products, and the cessation and/or reduction in the use of their synthetic counterparts [1]. Therefore, the growing demand and wide applicability of aroma compounds pose new challenges and require new trends in production. Conventional production procedures, such as hydrodistillation and extraction with organic solvents that are harmful to human health and the environment, cannot ensure the rational use of natural raw materials and a high-quality and safe product in adequate quantities in an environmentally benign production process. Additionally, aroma compounds are characterized by instabilities, easy reactivities, volatilities, and tendencies to degradation; therefore, their application is difficult [3]. Apart from organoleptic alterations and impairment of the compounds' biological properties, is it also possible to compromise the consumers' well-being, for example, by manifesting skin-sensitizing effects [3]. In order to prevent the deterioration of the quality and biological properties, it is necessary to ensure the stability of aroma compounds during the process of isolation, processing, and storage.
Furthermore, addressing the mentioned challenges is in accordance with the sustainable development agenda [4], as well as with the objectives of green chemistry, as it involves the creation of strategies that include the efficient use of natural renewable resources, attainment of long-term stable aroma compounds, elimination of the use of harmful solvents, waste reduction, and greater process efficiency.
For aroma preservation, the most industrially relevant and represented technique is spray drying. However, this technology implies high process temperatures, which can be a huge disadvantage for thermosensitive components. Also, numerous techniques have been developed for the preservation of aroma volatile compounds such as fluidized bed coating, freeze-drying, spray cooling and spray chilling, extrusion, co-crystallization, simple or complex coacervation, inclusion complexation, ionic gelation, the formation of liposomes, etc. [5,6]. In addition to ensuring high quality and preventing the loss of volatile compounds during obtaining and processing the product, a very important aspect is the simplicity and economy of the process. All preservation techniques represent additional processes after aroma attainment and require additional equipment. Therefore, the development of simple and economic procedures for obtaining and preserving aroma volatile compounds is necessary.
Supercritical CO 2 is a green solvent that has been introduced into industrial processes as a sustainable alternative to conventional hazardous organic solvents. Supercritical CO 2 is suitable for the efficient recovery of lipophilic aroma components from natural materials and enables obtaining extracts of a good safety, quality, and purity profile. When CO 2 is in a supercritical state (critical temperature 31.06 • C and critical pressure 7.39 MPa), by changing the process parameters, it is possible to adjust the thermophysical properties of the solvent and the characteristics of the product. In addition to being safe for work, cheap, and available, its advantage is that it is easily removed from the product and reused while the product does not need to be purified. Moreover, this results in reduced energy consumption, lower costs, and increased sustainability of the process. Additionally, degradation is avoided because the process is carried out at low temperatures [7,8].
The new generation of alternative solvents, natural deep eutectic solvents (NADES), comprise mixtures of components of natural origins that form intermolecular interactions when mixed in an appropriate molar ratio. The bonds they form are mainly through hydrogen bonding, which results in a significant depression of the melting point of the mixture and the formation of a liquid mixture. What distinguishes them is easy preparation with 100% atom economy, thermal and chemical stability, low-cost, non-flammability, biodegradability, and absence of or low toxicity. NADES mixtures are tailor-made solvents; hence, it is possible to tune their physico-chemical properties, making them versatile solvents in different fields of applications [9,10]. There are literature data that indicate the stabilizing properties of NADES for different types of components [11]. However, the potential of NADES for stabilizing aroma volatile organic compounds has not been investigated previously.
Therefore, the goal of this study was to establish a novel green approach of using environmentally friendly solvents to enable the attainment of aroma volatile compounds and maintain their stability at room temperature. Satureja montana L. (Lamiaceae) was selected for research, as it represents a rich source of volatile aroma compounds such as carvacrol, thymol, and thymoquinone, which are highly important and have applications in various fields. Due to its strong affinity for volatile compounds, supercritical CO 2 was used for the isolation of volatile compounds and obtaining clean and safe extracts, which were further stabilized using NADES. The stabilizing capacity of NADES mixtures was evaluated by monitoring the chromatographic changes in the volatile profile of the CO 2 extracts dispersed in different NADES during storage for 6 months. For data analysis and the assessment of the stability of the samples during time, the Wilcoxon matched pair test and different chemometric methods (hierarchical cluster analysis, principal component analysis, and a ranking approach-sum of ranking differences) were used with the aim of establishing adequacy of different NADES for the stabilization of S. montana volatile aroma compounds.

Materials and Chemicals
Satureja montana L. aerial parts were obtained from the Institute for Medicinal Plants Research "Dr. Josif Pancic", Belgrade, Serbia. The mean particle size (0.32 ± 0.05 mm) of the material was determined using vibration sieve sets (CISA, Cedaceria, Barcelona, Spain). Moisture content of the plant material (9.26%) was determined using a moisture analyzer DAB (Kern, Ballingen, Germany).

Preparation and Characterization of NADES
NADES mixtures were prepared by mixing components in an adequate molar ratio and then heating (40 • C) and stirring until a clear liquid was formed. In this work, different NADES based on betaine and glycerol were used. The NADES that were applied, abbreviation, and molar ratio are listed in Table 1.

Determination of Viscosity
The viscosity of the different NADES (Bet/Sor/W, Bet/Suc/Pro/W, Gly/Glu/Sorb/W, Pro/Gly/Sorb/W) was determined using a MCR102 Modular Compact Rheometer (Anton Parr) fitted with a parallel plate geometry with 50 mm of diameter (PP50, Anton Parr) and 1 mm of gap. Measurements were conducted in triplicates. Viscosity of Gly/Glu, Bet/EG, and Bet/Gly/Suc/W was determined previously [12,13].

Preparation of Samples
Supercritical Carbon Dioxide Extraction and Dispersion in NADES Supercritical CO 2 extraction was carried out using a high-pressure extraction system (HPEP, NOVA-Swiss, Effretikon, Switzerland) with the following main specifications: gas cylinder with CO 2 , diaphragm-type compressor (with a pressure range up to 1000 bar), extractor vessel with heating jacket (internal volume 200 mL, maximum operating pressure 700 bar), separator with cooling jacket (internal volume 200 mL and maximum operating pressure 250 bar), pressure control valve, temperature regulation system, and regulation valves. The extractions were performed under the following conditions: pressure 300 bar, temperature 40 • C, CO 2 flow 0.194 kg/h, and extraction time 4.5 h. The parameters of extraction were selected based on the previously published work [14]. Obtained extracts were placed in glass bottles and stored at 4 • C prior to further analysis. Each extraction was performed in triplicate.
After 4.5 h of extraction, to separate the CO 2 and extract, depressurization was performed, and the extract was directly dispersed into the prepared NADES mixture. The mixture of CO 2 extract and NADES was homogenized by mixing in a vortex for 1 min. The extract: NADES ratio was 0.05 g ± 15% extract: 1 mL NADES. This ratio was chosen because it enabled easy homogenization of the mixture. The control and CO 2 extracts that dispersed into NADES (CO 2 -NADES systems) were analyzed at the beginning (month 0), after 3 months, and after 6 months of storage in transparent containers at room temperature in the dark.

Headspace Solid-Phase Microextraction (HS-SPME)
Volatiles in the headspace were extracted using a solid phase microextraction (SPME) fiber installed on the PAL Auto Sampler System (PAL RSI 85, CTC Analytics AG, Schlieren, Switzerland). The fiber was covered with a layer of divinylbenzene/carbon wide range/ polydimethylsiloxane (DVB/Carbon WR/PDMS) and was purchased from Agilent Technologies (Palo Alto, Santa Clara, CA, USA). Before the analysis, the fiber was conditioned according to the manufacturer's instructions. The sample (1.5 g) was placed in a 15 mL glass vial and hermetically sealed with PTFE/silicone septum. The sample was equilibrated at 60 • C for 15 min, and extraction was performed for 45 min, followed by 6 min thermal desorption through the 250 • C inlet and directly onto the GC column.

Gas Chromatography with Mass Spectrometry Analysis (GC-MS)
An Agilent 8890 gas chromatograph (Agilent Technologies, Palo Alto, CA, USA) coupled to a mass spectrometer (series 5977E, Agilent Technologies, Palo Alto, CA, USA) was used for the analysis. The analyzed components were separated in a HP-5MS capillary column (30 m × 0.25 mm, 0.25 µm, Agilent Technologies, Palo Alto, CA, USA). The GC conditions were described previously [15].
Qualitative identifications of the compounds were performed using the Wiley 9 (Wiley, New York, NY, USA) and NIST 17 (National Institute of Standards and Technology, Gaithersburg, MD, USA) mass spectral libraries, as well as the literature data of retention indexes calculated with C 9 -C 25 alkanes. The HS-SPME/GC-MS analysis of each sample was performed in three replicates, and the results (relative percentages) were expressed as mean data.

Chemometric Methods
The chemometric analysis was performed on natural data since all the variables were on the same scale. HCA was carried out using the NCSS 2023 program [16]. The clustering of the samples and variables was presented on clustered heat maps (double dendrogram) to clearly show based on which variables the samples were grouped in the clusters. The clustering methods were Ward's variance, and the distance method was based on Euclidean distances. PCA was performed by using Statistica v. 14.0.0.15 program [17]. The analysis was based on correlations. The relevant number of PCs was estimated based on Eigenvalues higher than 1. The results of the PCA were presented in the form of scores and loadings plots that represented the distribution of the samples in the space of the considered variables. SRD analysis was performed by using a program created in Microsoft Excel software, created by Héberger and Kollár-Hunek [18,19]. The method was based on the calculation of the sum of ranking differences between the objects (samples, models, etc.) and a reference ranking. In the present study, the reference ranking was the row average values (consensus ranking). The SRD modeling was validated by Comparison of Ranks by Random Numbers (CRRN) approach and 7-fold cross-validation procedure [19].

Results and Discussion
Instability, volatility, reactivity, and inclination towards degradation are responsible for significant changes in aroma volatile components, which altogether has significant economic and health consequences. These deteriorations are particularly expressed at room temperature [3,20]. Therefore, storage at room temperature was chosen for stability testing. By ensuring the stability of aroma volatile compounds at room temperature, the application of these components is facilitated, the storage and transportation processes of such products are simplified, and health risks are reduced.
Easily available and low-cost constituents were selected for the design of NADES: glucose, sucrose, glycerol, sorbitol, proline, betaine, and ethylene glycol. The viscosity of NADES could represent an obstacle to applying NADES; however, it could potentially have a positive effect on the stabilizing capacity of NADES because it slows down the movement of molecules [21]. Therefore, NADES with wide viscosity ranges were used in this work (Table 1).
At the beginning (month 0), the headspace composition of S. montana extracts dispersed in different NADES (CO 2 -NADES systems) was established, followed by monitoring of the changes in the headspace profiles of the CO 2 -NADES systems after 3 and 6 months of storage at room temperature. CO 2 extract was used as a control, and its volatile profile was monitored in the same way.
The S. montana samples were found to be complex mixtures of components belonging to different chemical groups, such as hydrocarbons, alcohols, aldehydes, esters, ethers, ketones, phenols, and oxides. A total of 69 components were identified in the samples, which made up 83.06-98.97% of the headspace profile of the samples (Tables 2-4). Oxygenated terpenes were predominant in the samples and accounted for 53.99-86.05% of the volatile profile, followed by sesquiterpene hydrocarbons 3.29-28.03%, monoterpene hydrocarbons 1.01-9.38%, oxygenated sesquiterpenes 0.17-13.27%, and non-terpenes 0.40-3.97%. Tables 2-4 show the headspace composition of the control and CO 2 -NADES systems at the start (month 0), after 3 months, and after 6 months of storage.
Among the most dominant components in the samples were the isomeric monoterpene phenols carvacrol and thymol (6.38 ± 0.97-31.44 ± 3.21% and 5.65 ± 0.36-23.48 ± 1.77%, respectively). These two components are of exceptional commercial importance and have a wide spectrum of applications in other fields such as pharmacy, pharmacy, medicine, dentistry, cosmetics, food, veterinary, and agrochemistry [22]. According to the literature, S. montana lipophilic extracts and essential oils mainly belong to phenolic chemotypes with dominant carvacrol and/or thymol. Therefore, the distribution of components established in this work was in accordance with the literature [23,24].
Monitoring of the headspace profile (during 6 months) revealed changes in the chemical profile of the samples as a function of storage time. In the Bet/Gly/Suc/W sample, a constant decrease in thymol (from 20.84 ± 2.68 to 14.02 ± 3.14%) and carvacrol (from 25.56 ± 4.22 to 16.70 ± 2.50%) was observed during storage. Moreover, in the Bet/EG sample, thymol and carvacrol remained almost unchanged in percentage after the initial increase at 3 months, even after 6 months ( Figure 1). In other samples, carvacrol had an increasing trend (after 3 months) and then a decrease after 6 months. Thymol in the control and Gly/Glu exhibited the same trend-an increase (after 3 months) followed by a decrease. While, in the other samples (Bet/Suc/Pro/W, Bet/Sor/W, Pro/Gly/Sor/W, and Gly/Glu/Sor/W), thymol had small oscillations in abundance after 3 months, followed by a decrease (after 6 months), possibly indicating the preserved stability of thymol in the first three months in these samples.  The initial presence of carvacrol and thymol in the Bet/EG sample was significantly lower (6.38 and 5.65 ± 0.36%, respectively), compared to the remaining samples (thymol 16.25 ± 2.57-23.48 ± 1.77% and carvacrol 19.96 ± 3.04-29.37 ± 3.02%). Also, the distribution of other less abundant components in this sample was different from the remaining CO 2 -NADES systems and control. This difference could be explained by the different batches and non-uniformity of the material that was the commercial sample. Therefore, the monitoring of all samples from the beginning (month 0) to 6 months was of great importance.
Another important constituent that bettered the quality of the S. montana extract was thymoquinone. This component is classified as a promising candidate for drug development because it exhibits numerous pharmacological activities, high therapeutic indexes, favorable pharmacokinetics, and good safety and toxicity profiles [25]. Determined thymoquinone oscillations in the samples over time could potentially be a consequence of changes in the abundance of thymol and carvacrol because thymoquinone can be formed by subsequent oxidation of thymol and carvacrol. However, it was not possible to clearly establish these processes because most likely other components were included in the sequential oxidative reactions. For example, the changes in the presence of thymol, carvacrol, and thymoquinone could also be caused by changes in p-cymene, the precursor of thymol and carvarol. Furthermore, changes in p-cymene content could be preceded by other changes, such as the cyclization of myrcene to α-and γ-terpinenes and limonene, which may be further converted to p-cymene by rearrangement, hydrogenation, or dehydrogenation [26].  p-Cymene was dominant (0.84 ± 0.17-8.65 ± 2.58%) in the group of monoterpene hydrocarbons. In addition to having significant pharmacological properties [27], its presence in samples is essential because it is considered responsible for the more effective manifestation of the antimicrobial action of carvacrol. Namely, due to its hydrophobicity, it has a high preference for cytoplasmic membranes, allowing carvacrol to easily enter into the cell [28]. p-Cymene as a monoterpene hydrocarbon is more susceptible to changes compared to its oxidation derivatives. In the control, the stability of p-cymene was low compared to the CO2-NADES samples; so, its abundance in the control was reduced more than twice after 3 months, which coincided with the initial growth of carvacrol and thymol in the control. In the Bet/Gly/Suc/W sample, after 3 months, an increase in p-cymene of 2.8 times was recorded and, after that, a decrease to a percentage close to the initial one. In the other samples, p-cymene had significantly less oscillation over time, so its content was close to the initial one after 3 months. The highest abundance of p-cymene was in Bet-EG, where its percentage after 6 months (6.38%) was close to that at the start (7%).
Furthermore, a reduction in carvacrol methylether was observed in the samples, which may have been due to its degradation to carvacrol and p-cymene [29]. Linalool increased in all samples after 3 months, and it increased even after 6 months in most of them, with the exception of Bet/Gly/Suc/W and Bet/Sor/W, where it stagnated after the initial growth. Growth of linalool has been reported during storage of lemon-flavored hard tea at room temperature [30]. According to He et al. [30], the increase in linalool can occur due to autoxidation of myrcene [31], which, in the samples, was detected only at the beginning and also in Bet/EG in traces after 3 months. The increase in a minor compound αterpineol, which can be formed by oxidation of limonene, was also detected. The changes in the monoterpene alcohols terpinene-4-ol and borneol could also be attributed to the degradation reaction of hydrocarbons. p-Cymene was dominant (0.84 ± 0.17-8.65 ± 2.58%) in the group of monoterpene hydrocarbons. In addition to having significant pharmacological properties [27], its presence in samples is essential because it is considered responsible for the more effective manifestation of the antimicrobial action of carvacrol. Namely, due to its hydrophobicity, it has a high preference for cytoplasmic membranes, allowing carvacrol to easily enter into the cell [28]. p-Cymene as a monoterpene hydrocarbon is more susceptible to changes compared to its oxidation derivatives. In the control, the stability of p-cymene was low compared to the CO 2 -NADES samples; so, its abundance in the control was reduced more than twice after 3 months, which coincided with the initial growth of carvacrol and thymol in the control. In the Bet/Gly/Suc/W sample, after 3 months, an increase in p-cymene of 2.8 times was recorded and, after that, a decrease to a percentage close to the initial one. In the other samples, p-cymene had significantly less oscillation over time, so its content was close to the initial one after 3 months. The highest abundance of p-cymene was in Bet-EG, where its percentage after 6 months (6.38%) was close to that at the start (7%).
Furthermore, a reduction in carvacrol methylether was observed in the samples, which may have been due to its degradation to carvacrol and p-cymene [29]. Linalool increased in all samples after 3 months, and it increased even after 6 months in most of them, with the exception of Bet/Gly/Suc/W and Bet/Sor/W, where it stagnated after the initial growth. Growth of linalool has been reported during storage of lemon-flavored hard tea at room temperature [30]. According to He et al. [30], the increase in linalool can occur due to autoxidation of myrcene [31], which, in the samples, was detected only at the beginning and also in Bet/EG in traces after 3 months. The increase in a minor compound α-terpineol, which can be formed by oxidation of limonene, was also detected. The changes in the monoterpene alcohols terpinene-4-ol and borneol could also be attributed to the degradation reaction of hydrocarbons.
The group of oxygenated sesquiterpenes showed an increase in the samples, corresponding to a decrease in sesquiterpene hydrocarbons over time. In the group of oxygenated sesquiterpenes, a total of three components were identified, with an increasing trend over time in all samples compared to the start, which was caused by the transformations of sesquiterpene hydrocarbons. The appearance of caryophyllene oxide after 3 months and an increase up to 6 months could be attributed to the oxidation of β-caryophyllene, which decreased over time. An increase in caryophyllene oxide was reported in the literature as one of the characteristic indicators of the aging of essential oil and lipophilic products [3]. The identification of viridiflorol in the samples after 6 months was in accordance with the decrease and disappearance of the corresponding sesquiterpene aromadendrene hydrocarbon ledene (viridiflorene). Also, apart from viridiflorol, another aromadendron derivative, alcohol spathulenol, was identified.
In the control, a decrease in monoterpene hydrocarbons was observed approximately 3 times after 3 months. This trend in hydrocarbons in the control was in accordance with the literature's results, where a decreasing trend in low boiling hydrocarbons at room temperature storage was reported. Most likely, evaporation, oxidation, and other unwanted changes in the product constituents during the storage period were responsible for such changes [20]. On the other hand, in the CO 2 -NADES systems, the changes in the presence of hydrocarbons were significantly milder, mostly with a decreasing trend, potentially indicating that NADES improved the stability of these components compared to the control. The only CO 2 -NADES system in which a more pronounced oscillation (increase) was determined in the presence of this group of components was Bet/Gly/Suc/W.
Unstable constituents, such as hydrocarbons, are susceptible to chemical transformations due to their interactions with other components; hence, their stabilization is important for the stability of the entire product. As a result of the oxidation of the components, changes in the smell, taste, and texture of the product may occur. In addition, oxidation and other degradation reactions can lead to the formation of low molecular weight volatile components that exhibit undesirable effects on human health and lead to undesirable odors and flavors [32,33].
In the group of non-terpenes, several components were identified that were not initially present in the samples. The largest number of newly formed components, as well as the largest increase in the presence of non-terpene compounds, was observed in the control. In control after 6 months E,E-hepta-2,4-dienal was identified. This is a polyunsaturated aldehyde that can be developed due to a series of enzymatic transformations of mainly free polyunsaturated fatty acids. It is usually detected in stored oil samples with reduced stability and points to the oxidative changes related with lipid food spoilage and leads to the rancid defect [34]. Also, after 6 months, γ-butyrolactone, a decomposition product of tetrahydrofuran, was identified in the control, which was found to be acutely toxic at low concentrations in larval zebrafish [35]. Acetic acid was detected only in the control samples and over time its presence increased (up to 1.45%). This growth of acetic acid could indicate the reduced stability of the control. Huang et al. [36] stated that acetic acid can appear as a by-product of furan formation. Furthermore, the presence of acetic acid was reported as a consequence of undesirable processes during storage of honey and was accompanied by the development of an unpleasant and sour taste [37]. 4-Methylbenzaldehyde, the benzaldehyde derivative, was detected in all samples after 6 months, except for Gly-Glu and Bet/EG. Benzaldehydes are used as food preservatives; however, it has recently been established that prolonged exposure to benzaldehydes is associated with cancer and disorders of the liver, kidney, nervous, and reproductive systems [38].

Hierarchical Cluster Analysis and Principal Component Analysis
Statistical and chemometric methods were applied for further analysis and interpretation of the obtained results, which can be a useful tool in the analysis of large data sets [39][40][41]. In this work, the pattern recognition methods HCA and PCA, a ranking approach SRD, and Wilcoxon matched pair test were applied.
The HCA resulted in two double dendrograms. The dendrogram that represents the clustering of the control and CO 2 -NADES systems based on the percentage of each detected compound is given in Figure 2. The presented vertical dendrogram indicates the following findings: (1) the most abundant compounds in almost all extracts were carvacrol and thymol (placed in the purple vertical cluster), except in the Bet/EG extracts, in which thymoquinone was the most abundant compound; (2) the significant presence of borneol, terpinene-4-ol, caryophyllene oxide, linalool, carvacrol methyl ether, trans-caryophyllene, and p-cymene was observable in the majority of the extracts (placed in the separate vertical sub-cluster of the red-colored cluster); (3) the majority of the compounds were minor compounds, located in another sub-cluster of the vertical red cluster. The HCA resulted in two double dendrograms. The dendrogram that represents the clustering of the control and CO2-NADES systems based on the percentage of each detected compound is given in Figure 2. The presented vertical dendrogram indicates the following findings: (1) the most abundant compounds in almost all extracts were carvacrol and thymol (placed in the purple vertical cluster), except in the Bet/EG extracts, in which thymoquinone was the most abundant compound; (2) the significant presence of borneol, terpinene-4-ol, caryophyllene oxide, linalool, carvacrol methyl ether, trans-caryophyllene, and p-cymene was observable in the majority of the extracts (placed in the separate vertical sub-cluster of the red-colored cluster); (3) the majority of the compounds were minor compounds, located in another sub-cluster of the vertical red cluster. The horizontal dendrogram in Figure 2 shows the grouping of the CO2-NADES systems and control samples. It indicates that all three control samples were placed in different clusters and sub-clusters, suggesting that they were potentially different and, therefore, not very stable over time. Bet/EG samples were placed in a separate cluster (green The horizontal dendrogram in Figure 2 shows the grouping of the CO 2 -NADES systems and control samples. It indicates that all three control samples were placed in different clusters and sub-clusters, suggesting that they were potentially different and, therefore, not very stable over time. Bet/EG samples were placed in a separate cluster (green horizontal cluster), which indicated that they were distinctive from the other CO 2 -NADES systems, as previously explained. The main reason for their separation into different clusters was the content of thymoquinone, which was higher than in other CO 2 -NADES systems, and the significantly lower abundance of thymol and carvacrol in comparison with other extracts. The distribution of other NADES in the purple and red horizontal clusters showed that some samples monitored at the start and after 3 months were placed in the same cluster (e.g., Gly/Glu/Sor/W0 and Gly/Glu/Sor/W3, as well as Pro/Gly/Sor/W0 and Pro/Gly/Sor/W3, which belong to the red horizontal cluster), as well as some samples monitored after 3 and 6 months (Bet/Gly/Suc/W3 and Bet/Gly/Suc/W6). The main observations were scattering of the control samples in the dendrogram and significant separation of Bet/EG samples from other samples but in the same cluster, indicating that the Bet/EG samples had not suffered significant changes over time.
Another aspect of HCA of CO 2 -NADES systems and control is presented in Figure 3, where the clustered heat map shows the clustering of the samples in regard to the groups of compounds (oxygenated monoterpenes (OM), monoterpene hydrocarbons (MH), nonterpenes (NT), oxygenated sesquiterpenes (OS), and sesquiterpene hydrocarbons (SH)). The PCA modeling, based on the chemical composition in terms of the detected group of compounds, resulted in a model that took into account 76.5% of total variance based on two PCs described by Eigenvalues (Ei) higher than 1: PC1 (Ei = 2.51) explained 50.29%, and PC2 (Ei = 1.31) explained 26.20% of the total variance. In the score plot, presented in Figure 4a, it is evident that there was a clear separation (red line) between the samples monitored at the start of the experiment and the samples monitored after 3 and 6 months. The distribution of the samples along the PC1 axis was mostly influenced by SH, MH, and OM content, whereas the OS content had a significant influence on the PC2 axis (the score plot, Figure 4b). The content of NT had a moderate influence on the distribution of the compounds along both the PC1 and PC2 axes. The PCA results were in agreement with the results obtained by HCA presented in Figure 3. The main observations arising from the score graph indicated the distinction of the Bet/EG samples and significant differences between the control samples (significant distance of the sample Cont0 from the other control samples is evident). Additionally, the score plot indicated a noticeable change in the content of certain groups of compounds (SH, MH, and OM) from 0 to The vertical dendrogram indicates the significant separation of the OM compounds considering their highest contents (it can be considered as an outlier since it is placed outside the main cluster), whereas the other groups of compounds (MH, NT, OS, and SH) belonged to the same cluster considering their abundance, which was significantly lower in comparison to the OM compounds. The clustering of the CO 2 -NADES and control samples is presented in the horizontal dendrogram, which indicates the scattering of the control samples in different clusters and sub-clusters. Moreover, it also demonstrates the grouping of Bet/EG samples in a similar way as in the previous analysis, with the difference that they were placed in the same cluster (green horizontal cluster) with the CO 2 -NADES systems monitored at the start of the experiment (Gly/Glu0). The remaining CO 2 -NADES systems were placed in a separate cluster (red and purple horizontal cluster).
The PCA modeling, based on the chemical composition in terms of the detected group of compounds, resulted in a model that took into account 76.5% of total variance based on two PCs described by Eigenvalues (Ei) higher than 1: PC1 (Ei = 2.51) explained 50.29%, and PC2 (Ei = 1.31) explained 26.20% of the total variance. In the score plot, presented in Figure 4a, it is evident that there was a clear separation (red line) between the samples monitored at the start of the experiment and the samples monitored after 3 and 6 months. The distribution of the samples along the PC1 axis was mostly influenced by SH, MH, and OM content, whereas the OS content had a significant influence on the PC2 axis (the score plot, Figure 4b). The content of NT had a moderate influence on the distribution of the compounds along both the PC1 and PC2 axes. The PCA results were in agreement with the results obtained by HCA presented in Figure 3. The main observations arising from the score graph indicated the distinction of the Bet/EG samples and significant differences between the control samples (significant distance of the sample Cont0 from the other control samples is evident). Additionally, the score plot indicated a noticeable change in the content of certain groups of compounds (SH, MH, and OM) from 0 to 3 months for certain samples (e.g., control samples, Bet/EG, Gly/Glu, Bet/Suc/Pro/W, Bet/Sor/W).

The Ranking of the Samples Based on the SRD Approach
The ranking of the samples in terms of the order of observations (abundance of each detected compound) was carried out based on row average ranking as a golden standard (consensus ranking). The obtained results are presented in Figure 5. Based on the results, it can be concluded that there was no particular grouping of the samples regarding their SRD values. However, it can be noticed that the sample Gly/Glu/Sor/W3 was closest to the reference ranking, and the control samples Cont0 and Cont6 were placed furthest from the reference ranking, whereas the Cont3 sample was placed among the other CO2-NADES samples. This implied that the control samples Cont0 and Cont6 had a significantly different order of the compounds than the reference ranking and CO2-NADES systems. This could also imply that the variations in the change in composition in the control samples were more pronounced than in the case of the CO2-NADES samples. All the samples were significantly distant from the Gaussian curve, indicating that the results were not random in nature.

The Ranking of the Samples Based on the SRD Approach
The ranking of the samples in terms of the order of observations (abundance of each detected compound) was carried out based on row average ranking as a golden standard (consensus ranking). The obtained results are presented in Figure 5. Based on the results, it can be concluded that there was no particular grouping of the samples regarding their SRD values. However, it can be noticed that the sample Gly/Glu/Sor/W3 was closest to the reference ranking, and the control samples Cont0 and Cont6 were placed furthest from the reference ranking, whereas the Cont3 sample was placed among the other CO 2 -NADES samples. This implied that the control samples Cont0 and Cont6 had a significantly different order of the compounds than the reference ranking and CO 2 -NADES systems. This could also imply that the variations in the change in composition in the control samples were more pronounced than in the case of the CO 2 -NADES samples. All the samples were significantly distant from the Gaussian curve, indicating that the results were not random in nature. The right y-axis presents the relative frequencies for the Gauss curve (XX1 = 5% limit, Med = Median, XX19 = 95% limit).
The results of 7-fold cross-validation of the SRD procedure are presented in Figure 6 in the form of a Box-Whisker plot. The results are based on normalized SRD values (SRDnorm) obtained in each repeated ranking procedure. The CO2-NADES and control samples were ordered in the same way as in the SRD graph ( Figure 5). For each sample, there were seven SRD values obtained for each repeated SRD analysis. The median, 25%, and 75% of the range, minimum, and maximum The right y-axis presents the relative frequencies for the Gauss curve (XX1 = 5% limit, Med = Median, XX19 = 95% limit).
The results of 7-fold cross-validation of the SRD procedure are presented in Figure 6 in the form of a Box-Whisker plot. The results are based on normalized SRD values (SRDnorm) obtained in each repeated ranking procedure. The right y-axis presents the relative frequencies for the Gauss curve (XX1 = 5% limit, Med = Median, XX19 = 95% limit).
The results of 7-fold cross-validation of the SRD procedure are presented in Figure 6 in the form of a Box-Whisker plot. The results are based on normalized SRD values (SRDnorm) obtained in each repeated ranking procedure. The CO2-NADES and control samples were ordered in the same way as in the SRD graph ( Figure 5). For each sample, there were seven SRD values obtained for each repeated SRD analysis. The median, 25%, and 75% of the range, minimum, and maximum The CO 2 -NADES and control samples were ordered in the same way as in the SRD graph ( Figure 5). For each sample, there were seven SRD values obtained for each repeated SRD analysis. The median, 25%, and 75% of the range, minimum, and maximum values are presented as a Box-Whisker plot. The graph points to the consistency in the SRD analysis of the samples. An extreme SRDnorm value is observable for the sample Gly/Glu/Sor/W3; however, it can be considered insignificant. The vertical red dashed line shows the significant difference between the Cont0 and Cont6 samples in comparison to the other samples based on the Wilcoxon Matched Pairs Test at the 0.05 confidence level. This is a confirmation of the previously mentioned assumption that Cont0 and Cont6 possessed a significantly different order of the compounds content than the reference ranking and CO 2 -NADES samples.

Wilcoxon Matched Pairs Test of the NADES Samples
Pairwise comparison of samples observed at the beginning of the experiment and after 3 and 6 months was performed using the Wilcoxon Matched Pairs Test to test whether the mean values of the samples differed over time (Table 5). The significant composition change between 3 and 6 months was noticeable for the control, Bet/Suc/Pro/W, Bet/Gly/Suc/W, Pro/Gly/Sor/W, and Gly/Glu/Sor/W samples. The Wilcoxon test implied the stability of the Gly/Glu, Bet/Sor/W and Bet/EG samples for 6 months. The improved stability of the components dispersed in individual NADES compared to the control after 6 months of storage at room temperature is evident. A possible explanation for the better stability of the CO 2 -NADES systems compared to the control is that NADES played the role of a physical protector of the components, which, consequently, ensured the protection of the components from oxidative degradation. Potentially, the high viscosity of Gly/Glu (5258.47 ± 322.86 mPa·s) and Bet/Gly/Suc/W (2246.77 ± 45.27 mPa·s) could be one of the factors that contributed to the stabilization. In a previous study, it had been suggested that stabilization in high-viscous NADES can contribute to the preservation due to the reduced movement of the molecules [21]. However, in our study, no clear relationship between viscosity and stability was noticed because Bet/EG (48.63 ± 0.42 mPa·s) with a very low viscosity proved to be one of the NADES with stabilizing capacity. Furthermore, the stabilizing capacity of NADES may be due to the development of intermolecular interactions with individual components of S. montana [11,21]. Namely, S. montana components can play the role of hydrogen bond donors and acceptors; hence, they can potentially form interactions with NADES components and thus slow down the movement of molecules, preventing oxidation and other degradation reactions and maintaining the overall stability. Also, for the further implementation of the established process on a large scale, an important aspect is the price of the used NADES (Table 1). Amino-acid-based NADES are evidently more expensive compared to sugar-based ones. Among the systems for which stability was established over time, the cheapest was Bet/Sor/W, followed by Bet/EG and Gly/Glu.
Conventional procedures for obtaining and stabilizing aroma volatile components include Soxhlet or hydrodistillation, followed by spray drying. Soxhlet extraction involves the use of hazardous organic solvents, as well as the possibility of their remaining in the final product. On the other hand, hydrodistillation is performed at high temperatures that can lead to the degradation of thermally labile compounds, and, as a result, the smell of hydrodistilled essential oil is often not authentic to the source. With the application of supercritical CO 2 at mild temperatures, the possibility of deterioration is reduced, there is no use or generation of toxic waste, and the obtained extract is clean and safe [42]. Furthermore, the spray drying technique is the most common method of preserving volatile components; however, it implies the use of carriers and additional processing of the sample, which, apart from additional costs, also suggests working at high temperatures and the possibility of losing thermo-sensitive components [5]. On the other hand, the preparation of NADES is simple, and dispersion in NADES avoids additional processing of samples at high temperature.
The resulting systems represent products with numerous potential applications. Due to the presence of components with important pharmacological properties, such as carvacrol, thymol, thymoquinone, and others, they can be used in obtaining pharmaceutical products. Also, the components of S. montana are highly relevant in the cosmetic industry as well as the perfume industry. Additionally, the properties of S. montana components, such as the improvement of sensory characteristics and preservative and antimicrobial properties, are also important for the food industry. Moreover, the obtained systems can play an important role in creating repellent products [43,44].

Conclusions
The approach investigated in this work represents the sustainable aroma stabilization process, which ensures the efficient attainment of aroma compounds and their preservation, consequently reducing costs, promoting a more rational use of natural resources, and reducing waste generation. Also, the goals of this work correspond to the objectives of green chemistry and the Sustainable Development Goals agenda, as they include environmental and economic sustainability, environmental protection, and human well-being.
By comparison, between the changes in the headspace profiles of the control and CO 2 -NADES systems, it is evident that the application of NADES advanced the stability of S. montana volatile compounds. Apart from the more pronounced changes with regard to oxidation and other transformational reactions, the formation of non-terpene components such as acetic acid, E,E-hepta-2,4-dienal, and γ-butyrolactone was observed in the control. Moreover, these components can significantly reduce the quality of the products. In CO 2 -NADES systems, there was no formation of these components, and the changes in the headspace profiles were less pronounced compared to the control. Bet/EG was observed as a system with the lowest oscillations over time. In addition, using the Wilcoxon test, it was established that the systems with the least changes were Bet/EG, Gly/Glu, and Bet/Sor/W.
It is of great importance to continue the study of aroma compounds from different types of materials and their interactions with different NADES because the stabilities of volatile components, essential oils, extracts, and other products containing volatile aroma compounds depend on the characteristics of those components and their tendencies to oxidation and degradation changes.