Weather Conditions Inﬂuence on Hyssop Essential Oil Quality

: This paper is a study of the chemical composition of Hyssopus ofﬁcinalis ssp. ofﬁcinalis grown during three years (2017–2019) at the Institute of Field and Vegetable Crops Novi Sad (Vojvo-dina Province, Serbia). Furthermore, comparisons with ISO standards during the years were also investigated, as well as a prediction model of retention indices of compounds from the essential oils. An essential oil obtained by hydrodistillation and analysed by GC-FID and GC-MS was isopinocamphone chemotype. The gathered information about the volatile compounds from H. ofﬁcinalis was used to classify the samples using the unrooted cluster tree. The correlation analysis was applied to investigate the similarity of different samples, according to GC-MS data. The quantitative structure– retention relationship (QSRR) was also employed to predict the retention indices of the identiﬁed compounds. A total of 74 experimentally obtained retention indices were used to build a prediction model. The coefﬁcient of determination for the training cycle was 0.910, indicating that this model could be used for the prediction of retention indices for H. ofﬁcinalis essential oil compounds.


Introduction
Hyssopus officinalis L. commonly known as hyssop, is a medicinal, culinary, essential oil-bearing, melliferous and ornamental plant [1,2]. This perennial plant with a strongbranching taproot grows as a subshrub, up to 70 cm high, with opposite shiny dark green lanceolate or oblong leaves and blue-white or violet flowers arranged in false spikes [3]. During the flowering stage, the aboveground parts (Hyssopi herba) taste slightly bitter and have a strong aromatic flavour, somewhat like a cross between sage and mint [4]. Hyssopus officinalis L. is mostly used for essential oil production. Its essential oil is used in the pharmaceutical and perfume industries and cosmetics as well as in aromatherapy [2][3][4][5].
It originates from the Mediterranean region; however, it is well adapted to the northern hemisphere [5]. This plant has two main subspecies in Serbia: ssp. aristatus [syn. ssp. pilifer] (in nature on limestones, dry and sunny slopes, and meadows) and ssp. officinalis (in plantations) [2,6].
The wild-grown subspecies contains linalool as a main compound of the oil, from (35.3 to 51.2)% [7]. This subspecies does not satisfy the requirements of ISO 9841 standard which does not recognise linalool as the main component of hyssop essential oil and does not have it listed or defined required concentration range. Ref. [8], as well as var. decumbens

Soil Characteristics
The soil used for plants growing in this experiment belongs to the chernozem calcareous gleyic type. Basic chemical properties of the soil were determined in samples of the top layer (0 to 30) cm. The soil samples were analysed at the Faculty of Agriculture, University of Novi Sad, in accordance with standardized methods adopted in Serbia. Soil pH value was determined potentiometrically in suspension with 1M KCl and suspension with distilled H2O (SRPS ISO 10390:2007) [20]. CaCO3 content was determined using the volumetric method (JUS ISO 10693:2005) [21]. The total nitrogen content was determined using the Semimikro-Kjeldahl method modified by Bremner [22]. Available phosphorous and potassium in soil were determined by extraction with ammonium lactate solution (AL) in accordance with the method established by Egner, Riehm and Domingo [23]. Results of agrochemical analyses are shown in Table 1.  Weather conditions in the third growing year (2018/2019): autumn (X-XI) was drier and hotter than usual, winter (XII-II) was hotter with the usual amount of precipitation. The beginning of spring was dry and hot with frost occurrence during March and April, as well as unstable weather with frequent large amounts of precipitation. Summer (VI-VIII) was hot (23.4 • C) and humid (174 mm). Weather conditions in the third growing year (2018/2019): autumn (X-XI) was drier and hotter than usual; winter (XII-II) was hotter with the usual amount of precipitation. The beginning of spring was dry and hot with frost occurrence during March and April, as well as unstable weather with frequent large amounts of precipitation. Summer (VI-VIII) was hot (23.4 • C) and humid (174 mm). As the last month of the vegetation period, September was categorized as favourable with hot, dry and sunny weather (important for next growing year).

Soil Characteristics
The soil used for plants growing in this experiment belongs to the chernozem calcareous gleyic type. Basic chemical properties of the soil were determined in samples of the top layer (0 to 30) cm. The soil samples were analysed at the Faculty of Agriculture, University of Novi Sad, in accordance with standardized methods adopted in Serbia. Soil pH value was determined potentiometrically in suspension with 1M KCl and suspension with distilled H 2 O (SRPS ISO 10390:2007) [20]. CaCO 3 content was determined using the volumetric method (JUS ISO 10693:2005) [21]. The total nitrogen content was determined using the Semimikro-Kjeldahl method modified by Bremner [22]. Available phosphorous and potassium in soil were determined by extraction with ammonium lactate solution (AL) in accordance with the method established by Egner, Riehm and Domingo [23]. Results of agrochemical analyses are shown in Table 1. GC-FID and GC-MS analyses were performed according to our previous work [24] with the exception of split ratio which was 10:1 in the present work. The compounds were identified based on a comparison of their linear retention index relative to C 8 -C 32 n-alkanes and mass spectra with MS data library (Adams4 and NIST17). The relative mass percentage concentrations were computed from GC-FID peak areas.

QSRR Analysis
The molecular structure dataset was presented in short ASCII files (with extension. smi, the simplified molecular-input line). The molecular data were collected from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/, accessed on 30 January 2021). The calculation of the specified molecular descriptors for each chemical compound [25] was performed using free molecular descriptor software PaDel-descriptor [26]. A huge amount of data was obtained for each chemical compound, and it was necessary to use a genetic algorithm (GA) [27,28], and the calculation was performed using Heuristic Lab (https://dev.heuristiclab.com/trac.fcgi/, accessed on 30 January 2021), to select the most relevant molecular descriptors for retention time prediction. In this study, GA was used to determine the most appropriate molecular descriptors to develop a reliable model for predicting the retention time of the compounds found in H. officinalis essential oil. The correlation between the descriptors was examined, and collinear descriptors were detected using factor analysis.

Artificial Neural Network (ANN)
A multi-layer perceptron model (MLP), which is explained by the input, hidden and output layers, was used in this paper, considering that it is well known as capable of approximating nonlinear functions [29]. Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm was used for ANN modelling. The input and output data were normalized to improve the behaviour of the ANN. The experimental database was randomly divided into training (60%), testing (20%) and validation parts (20%) for ANN modelling. ANN results, including the weight values, depending on the initial assumptions of parameters necessary for ANN construction and fitting [30,31]. A series of different topologies were used, in which the number of hidden neurons varied from 5 to 20 and the training process of the network was run 100,000 times with random initial values of weights and biases. The optimization process was performed based on validation error minimization. Statistical investigation of the data has been performed mainly by the Statistica 10 software [32].

Global Sensitivity Analysis
Yoon's interpretation method was used to determine the relative influence (RIN) of molecular descriptors on retention time [33]. This method was applied based on the weight coefficients of the developed ANN.

Cluster Analysis
The cluster analysis (CA) was employed to evaluate the variations of the EO's compounds content of H. officinalis samples gathered in different sites and/or taken from literature reports. The phylogenetic tree diagram for obtained samples was plotted using the R software 4.0.3 (64-bit version). The R package "ape" (Analysis of Phylogenetics and Evolution) was utilized for calculation ("unrooted" tree cluster analysis was plotted). The distance matrix was calculated using the Euclidean method, using the "complete" method.

Linear Regression Model
A linear regression model was developed to predict the H. officinalis essential oil active compound content according to temperature and precipitation data, and the appropriate regression coefficients were calculated.

Correlation Analysis
The correlation analysis was performed using the R software, to present the similarities between different samples graphically.

Statistical Analyses
The collected data were processed statistically using the software package STATISTICA 10.0 [32].
On isopinocamphone as the most abundant component (41.2% on average for three years), the influence of temperatures was negative (regression coefficient was −26.218), while the impact of precipitation was positive (1.336) on the accumulation. This effect explains the reduced accumulation in the 2nd year of the study (38.8%) when the mean daily temperatures were higher and the total rainfall during the growing year compared to the other years of the research. On the other hand, the influence of increased temperatures was positive on the accumulation of pinocamphone (regression coefficient was 20.860) and β-pinene (7.006), while the impact of precipitation was negative (−1.065 and −0.693, respectively).
According to the ISO 9841 Standard [8], the content of pinocamphone, isopinocamphone, β-caryophyllene, germacrene D, elemol and spathulenol satisfied requirements during all three years, while α-pinene, β-pinene, sabinene, limonene, β-bourbonene and allo-aromadendrene was slightly exceeded in some years. Myrtenyl methyl ether was not detected in this investigation (Table 3). Therefore, the high quality of H. officinalis essential oil indicates that there is interest in growing this species in plantations in Serbia.

Climatic Factors and Their Influence on Hyssop Essential Oil Chemical Composition
The chemical composition variations of H. officinalis oils can be related to many factors. Still, the most important is intraspecific diversity, as a consequence of adaptive evolution on diverse habitats (geographical and climatic) and technological influences (the stages of development and harvesting time, the part used, irrigation) [36,37]. However, our previous studies were dealing with chemical composition and biological activities [2,34,35]. It is known that biological activity depends on the chemical composition of essential oil [38,39]. During the three-year experimental period (2012-2014), the vegetation year significantly affected essential oil yield but not the essential oil composition [40]. Also, during three years of cultivation of hyssop in Finland, the content of pinocamphone varied between 18.0 and 19.7%, isopinocamphone between 29.2 and 32.6%, germacrene D between 6.7 and 10.2%, β-pinene between 5.2 and 6.6% and pinocarvone between 5.3 and 5.7% [41]. In addition, investigations conducted in Ukraine for three years showed significant variations in the content of pinocamphone (from 35.49 to 53.73%) and isopinocamphone (4.66 to 44.43)% depending on plant age [42].
A strong negative correlation was observed between the amount of precipitation and the accumulation level of volatile organic compounds in the flowering shoots: rain during the development of flowers seems to have a negative effect on the yield [40]. Trans-pinocamphone was predominant in the oil during the vegetative stage, however, its content decreased with plant growth, while the content of pinocamphone increased [43], gradually from the pre-blooming stage (47.94%) through the full-blooming (48.22%) to the post-blooming stage (51.42%) [11]. This phenomenon indicated that choosing a suitable harvesting stage to achieve the highest essential oil quality was crucial [11]. It is known that saturated bicyclic monoterpene ketones pinocamphone and isopinocamphone pathway rise from β-pinene via allyc alcohol trans-pinocarveol by subsequent oxidation and two stereochemical alternatives for reduction of the conjugated double bond [44].
This plant is a typical xerophyte and is well adapted to drought [11]. Furthermore, it is also cold and frost tolerant [45]. Draught stress increases the number of essential oil compounds (from 27 in non-stress plants to 42 in intense drought stress conditions). In addition, with water deficit, the content of cis-pinocamphone is decreased [46]. In semiarid conditions, with 300-400 mm annual rainfall, the content of isopinocamphone was 19.8% and pinocamphone was 17.2%, while in irrigated crops contents of both compounds increased to 23.6% and 18.5%, respectively [47]. However, severe drought citrulline could increase the yield of isopinocamphone, so it could be a novel metabolite agent used for improving qualitative parameters in hyssop cropping systems [5]. In addition, salt stress influences the relative content of components in essential oil. Pinocarvone content was significantly decreased with increasing NaCl levels, while isopinocamphone being the least affected [48]. Table 4 is shown 70 accessions of H. officinalis essential oil from literature, with 25 most abundant compounds (on average more than 0.5%). According to the unrooted cluster tree (Figure 2), it could be said that there are several chemotypes. Still, they mainly contain a mixture of pinocamphone + isopinocamphone + β-pinene in different proportions but also appear as a particular chemotype.      As β-pinene chemotype are classified two accessions (ssp. aristatus) with 19.3 and 24.7% of β-pinene [70]. As pinocamphone chemotype are classified accessions with (34.0 to 58.3)% of pinocamphone [35,41,42,[64][65][66]69]. However, this is much more than the requirements of ISO 9841 Standard (8.0 to 25.0)%.
The obtained data determine the position of the samples in the factor planes. The geometrically close points suggest the similitude of patterns represented by these points.
The obtained data determine the position of the samples in the factor planes. The geometrically close points suggest the similitude of patterns represented by these points.

Correlation Analysis
The correlation analysis was performed to investigate the similarities in active compounds content of the different Hyssopus officinalis samples , and the results were visually displayed in Figure 3. It can be noticed from the figure that the darker blue colour of the squares, which shows the two samples relation, presents a stronger correlation between these samples, i.e., the more pronounced similarity in the active compound content. On the other hand, the lighter tone suggests a certain dissimilarity between samples. Therefore, if the colour tone is slighter, the correlation is lower.

QSRR Analysis
The PaDel-descriptor was employed to evaluate the molecular descriptors. They present many aspects of the investigated compounds and have been successfully used in QSRR investigation. Prior to the GA calculation, the factor analysis was performed to eliminate the descriptors with equal or almost equal values of the examined molecular structures. Only one of the correlated descriptors remained in the GA calculation. GA was used to select the most appropriate set of molecular descriptors while choosing the most relevant set of descriptors was realized using the evolution simulation [80,81]. The number of genes was equal to the number of the molecular descriptors obtained in the

Correlation Analysis
The correlation analysis was performed to investigate the similarities in active compounds content of the different Hyssopus officinalis samples , and the results were visually displayed in Figure 3. It can be noticed from the figure that the darker blue colour of the squares, which shows the two samples relation, presents a stronger correlation between these samples, i.e., the more pronounced similarity in the active compound content. On the other hand, the lighter tone suggests a certain dissimilarity between samples. Therefore, if the colour tone is slighter, the correlation is lower.

QSRR Analysis
The PaDel-descriptor was employed to evaluate the molecular descriptors. They present many aspects of the investigated compounds and have been successfully used in QSRR investigation. Prior to the GA calculation, the factor analysis was performed to eliminate the descriptors with equal or almost equal values of the examined molecular structures. Only one of the correlated descriptors remained in the GA calculation. GA was used to select the most appropriate set of molecular descriptors while choosing the most relevant set of descriptors was realized using the evolution simulation [80,81]. The number of genes was equal to the number of the molecular descriptors obtained in the PaDel-descriptor. The population of the first generation in the GA calculation was selected randomly. The number of the genes was kept low during the GA calculation, to maintain a small subset of descriptors. The probability of generating zero for a gene was set at least 60%. The operators used in the simulation were: crossover (90% probability) and mutation (0.5%). A population size of 100 genes was chosen for GA, and evolution was allowed for over 50 generations. The evolution of the generations was stopped when 90% of the generations took the same fitness. As a result, the five most significant molecular descriptors selected by GA were: ATSC5e, AATSC1v, GATS3p, IC0 and MW.
Detailed explanations about the descriptors were found in the Handbook of Molecular Descriptors [82]. These descriptors encoded different aspects of the molecular structure and were applied to develop a QSRR model. Table 5 represents the correlation matrix among these descriptors. The calibration and predictive capability of a QSRR model should be tested through model validation. The most widely used squared correlation coefficient (r 2 ) can provide a reliable indication of the fitness of the model. Thus, it was employed to validate the calibration capability of a QSRR model.
To explore the nonlinear relationship between RIs and the selected descriptors, the ANN technique was used to build a retention time predictive model. The statistical results of the MLP 5-10-1 network are shown in Table 6. The quality of the model fit was tested in Table 4, in which the lower reduced chisquare (χ 2 ), mean bias error (MBE), root mean square error (RMSE), mean percentage error (MPE). Values show a better fit to the experimental results [83].
The predicted RIs are presented in Figure 4a to confirm the adequate prediction capabilities of the constructed ANN, by showing the relationship between the predicted and experimental retention values.
Obtained results reveal the reliability of the ANN models for predicting the RIs of compounds in H. officinal is essential oil determined by GC-MS.
In this section, the influence of the five most important input variables, identified using the genetic algorithm on RI was studied. According to Figure 4b, MW was the most influential parameter with an approximately relative importance of 67.11%, while the influence of GATS3p was −16.27%.

Conclusions
The results of the chemical composition analysis of Hyssopus officinalis ssp. officinalis grown during three growing years reviled isopinocamphone as the most abundant component, on which accumulation, the temperature had negative, while precipitation had a positive impact. The opposite effect was noticed on the content of pinocamphone and β-pinene. The chemical composition variations of H. officinalis oils can be associated with many factors; however, the most prominent is intraspecific diversity.
Obtained the excellent quality of H. officinalis essential oil indicated the possibility of producing this species in plantations in Serbia. Established the quantitative structure retention relationship model for predicting the retention times of chemical compounds in H. officinalis essential oil obtained by hydrodistillation and analysed by GC-MS using the genetic algorithm variable selection method and the artificial neural network model was shown as reliable and has great potential for practical use. The global sensitivity analysis-Yoon interpretation method identified molecular weight as the most influential parameter affecting the retention indices of chemical compounds in H. officinalis essential oil.  In this section, the influence of the five most important input variables, identified using the genetic algorithm on RI was studied. According to Figure 4b, MW was the most influential parameter with an approximately relative importance of 67.11%, while the influence of GATS3p was −16.27%.

Conclusions
The results of the chemical composition analysis of Hyssopus officinalis ssp. officinalis grown during three growing years reviled isopinocamphone as the most abundant component, on which accumulation, the temperature had negative, while precipitation had a positive impact. The opposite effect was noticed on the content of pinocamphone and β-pinene. The chemical composition variations of H. officinalis oils can be associated with many factors; however, the most prominent is intraspecific diversity.
Obtained the excellent quality of H. officinalis essential oil indicated the possibility of producing this species in plantations in Serbia. Established the quantitative structure retention relationship model for predicting the retention times of chemical compounds in H. officinalis essential oil obtained by hydrodistillation and analysed by GC-MS using the genetic algorithm variable selection method and the artificial neural network model was shown as reliable and has great potential for practical use. The global sensitivity analysis-Yoon interpretation method identified molecular weight as the most influential parameter affecting the retention indices of chemical compounds in H. officinalis essential oil.