Are environmental factors responsible for essential oil chemotype distribution of Balkan Juniperus communis var. saxatilis populations?

Abstract Juniperus communis var. saxatilis is a short prostrate evergreen conifer that grows on dry hills in the mountain areas. In the Balkans, this taxon has a fragmented distribution. Both leaves and cones are rich in essential oils (EOs), which play an important role in plant interaction with the environment. Even though leaves and berries represent a significant source for humans and wildlife, there is no research on the genetic and EO variability in the region. In this article, we analysed individuals from 18 populations from the entire Balkan Peninsula. Three bioclimatic regions were identified, based on precipitation and temperatures. Amplified fragment length polymorphism (AFLP) data revealed high genetic variability, which is mostly located within the populations (91.7%). In total, 156 compounds representing 99.0% of total oil were detected and identified. Three chemotypes (ɑ-pinene, sabinene and δ-3-carene) were detected in all samples. Although all chemotypes could be detected throughout the region, their distribution was correlated with environmental factors. The influence of environmental variables on EO composition and chemophenetics are discussed.


Introduction
Juniperus communis L. (Cupressaceae, section Juniperus) is a dioecious species with a wide, holarctic distribution. It grows on different substrates, at altitudes ranging from sea level to 2400 m.a.s.l. (Tutin et al. 1976;Adams 2011;Farjon and Filer 2013) across North America and Eurasia. Such a wide geographic and edaphic distribution is possible due to its very high genetic and phenetic plasticity. Another part of its success lies in fleshy seed cones . In short, the high genetic variability, coupled with the high phenotypic plasticity and rapid expansion of this taxon after the last glaciation, has attributed to taxonomic complexity and confusion within this taxon. In the past few decades, taxonomic relations have changed several times. Some authors separated this species into several closely related species (e.g. J. nana Willd., J. sibirica Burgsd., J. communis L., J. intermedia Schur., J. pygmaea K.Koch, etc.), while recent molecular data [both Random Amplified Polymorphic DNA (RAPDs) and chloroplast and nuclear sequences] have led to the treatment of all these taxa as varieties of a single species -J. communis L. (Adams 2011;Farjon and Filer 2013). Unlike typical J. communis variety, J. communis var. saxatilis Pall. (syn. J. nana Willd., J. sibirica Burgsd., J. montana (Aiton) Lindl. & Gordon, J. pygmaea K .Koch) (http://www.ipni.org 2021) is a short prostrate evergreen conifer that grows up to 50 cm. It grows on dry hills or mountain tracts, from 600 to 2400 m above sea level (Tutin et al. 1976;Adams 2011). In the Balkans, this variety has a fragmented distribution and can only be found at higher altitudes (Vidaković 1982;Adams 2011;Farjon and Filer 2013). The branches are short and thick, and the leaves are in whorls of three, ascending, lanceolate, 5-10 mm long, and abaxially keeled. Female plants bear fleshy seed cones (juniper berries or galbuli) usually 7-10 mm in diameter that are dark blue when ripe. Both leaves and cones are rich in essential oils (EOs) (Vidaković 1982;Adams 2011).
EOs are complex mixtures of volatile compounds produced by aromatic plants with many biological functions, both for plants and human use (Rajčević et al. 2013(Rajčević et al. , 2018(Rajčević et al. , 2020. These complex mixtures are under strict genetic control, thus showing their chemophenetic significance (Adams 2000;Rajčević et al. 2013Rajčević et al. , 2015Rajčević et al. , 2018Roma-Marzio et al. 2017;Dodoš et al. 2019). The latest investigation of the north-western populations of J. communis L.var. communis (Rajčević et al. 2020) has shown a separation of continental and coastal populations based on the EO composition, suggesting that environmental factors play a role in the distribution of EO chemotypes. A previous investigation of Corscian populations of this taxon by Ottavioli et al. (2009) revealed the existence of two different chemotypes -a rare sabinene chemotype, and a more common limonene chemotype, though there was no apparent pattern of distribution in such a small region.
Junipers and their berries are an important food source for birds and mammals, but also an irreplaceable crop in the US$12bn worth gin industry (Farris et al. 2017;Stricklan 2019;Stricklan et al. 2020;Gin -worldwide|Statista Market Forecast 2021). In Serbia, the cones are often added to a particular type of brandy ("rakija") called "Klekovača" and used to stimulate appetite (Lesjak et al. 2013). Both leaves and berries are also used in traditional medicine to treat respiratory diseases, cough, bronchitis and asthma (Vasilijević et al. 2018). Several studies indicated antioxidant, anti-inflammatory and antimicrobial activity (Judžentienė 2019;Darwish et al. 2020). All these bioactivities, but the aroma as well, are directly linked to the EO composition. For example, α-pinene is a highly bioavailable compound with rapid metabolism or redistribution, with a myriad of biological activities, e.g. antibacterial, antifungal, neuroprotective, gastroprotective, antitumor, that adds a woody aroma (Allenspach and Steuer 2021). On the other hand, sabinene possesses antifungal and anti-inflammatory activity and, unlike α-pinene, it imparts spicy odours used in the culinary and perfume industries (Cao et al. 2018). In addition, these compounds are still very hard to synthesise, so most of them are still being isolated from the various plant materials (Cao et al. 2018).
The previous investigation of the variability of this taxon, using the composition of cuticular n-alkanes, has shown a separation of continental from coastal populations (Rajčević et al. 2014) that was correlated with the climate, however only phytochemical markers were used in this study. Even though these results correspond to the ones obtained by Dodd and Poveda (2003), due to the complex natural history of the Balkan populations during the glaciations and interglaciations, without taking into account genetic variability, it is hard to fully attribute this divergence solely on the bioclimatic factors. Earlier analyses of a sister taxon, J. deltoides R.P. Adams, in the Balkans and Italy (Rajčević et al. 2015;Roma-Marzio et al. 2017) have also shown differentiation in EO composition of juniper depending on location, suggesting that the distribution of chemotypes is related to the climate, thus differentiating coastal and continental populations. With the exception of a few studies, there was never a systematic approach in the investigation of the EO composition and its geographic distribution, especially in relation to different pedological and environmental factors of this taxon. In addition, previous investigations of the genetic variability of this taxon included only RAPD molecular markers and a limited number of individuals (Adams 2000;Adams and Nguyen 2007). Amplified fragment length polymorphism (AFLP) was used in three different areas in Europe, but focusing mostly on J. communis var. communis, with no studies on var. saxatilis. These studies suggested high genetic diversity and low population differentiation (Merwe et al. 2000;Michalczyk et al. 2010;Vanden-Broeck et al. 2011). Given the changing climate that disproportionately affects high mountain areas where this variety grows, and considering that only a tiny amount of data is available on the EO composition of J. communis var. saxatilis throughout its vast range, as well as the importance of this taxon in industry and for local communities, the goal of this study was to analyse the variability of leaf EOs throughout the Balkan Peninsula, taking into account both genetic and environmental factors. The significance of the results of this study could be regarded from the ecologic and chemophenetic viewpoints. To the best of our knowledge, this combined approach using genetic markers (AFLP) and EO data was never used in studying this taxon.

Plant material
Twigs with needles of 376 individuals of J. communis var. saxatilis, from 18 natural populations were collected in (summer or autumn) from 2009 to 2015, along an altitudinal gradient from 1250 to 1950 m a.s.l. in the Balkan Peninsula. All samples were placed in individual plastic zip bags and put in a portable freezer on the location site. Within 24 h, the bags were transferred to the freezer and kept at a temperature of -18 °C to preserve EO composition until extraction. Plants from several locations were collected in different seasons to account for possible interseasonal variability. Plant material for DNA isolation was desiccated in the field using silica gel and stored at room temperature until DNA extraction. Population details, geographic information, terrain inclination, exposure, altitudes, collector names and voucher numbers are given in Table 1. The EO analysis included all 376 specimens, while 90 individuals were used for AFLP analysis. Ten individuals from the Orjen population were mixed in the field, so only one mixed sample of EO was obtained from this population. Taxonomic identification of the species was conducted by the authors in the field (NR and PM) (Barkworth and Adams 1993;Adams 2011). Voucher specimens were deposited in the BEOU (The herbarium of the University of Belgrade -Faculty of Biology).

Essential oil extraction
The needles (leaves) from the collected material were manually separated and ground prior to EO extraction. The homogenised material of each sample (ca. 5 g fresh weight) was subjected to 2-h simultaneous distillation and extraction in a Likens-Nickerson type apparatus using 5 mL of dichloromethane (Ch 2 CL 2 ) (Chaintreau 2001). Obtained extracts (0.5 mL) were stored in amber vials at 4 °C until further analyses.

Essential oil analysis
The EO composition was analysed using an Agilent 7890A apparatus equipped with a 5975C mass-selective detector, flame ionization detector (FID) and DB-5 MS fused-silica gel cap. The conditions and analysis procedures for GC-FID and GC/MS (Gas Chromatography -Mass Spectrometry) are described by Rajčević et al. (2020). In all experiments, the relative amounts of volatile components were expressed as percentages of the peak area of total ion chromatograms.
Values under 0.05% were not considered during compound identification. Library search and mass spectral deconvolution and extraction were performed using the software NIST AMDIS version 2.64.113.71, with the retention index (RI) calibration data analysis parameters set to "strong" level and 10% penalty for compounds without RI. The search was performed against our homemade library, containing 4972 spectra. The relative contents of the identified compounds were computed from the GC peak areas. Linear RI was calculated for all compounds using the following formula: LRI = 100 × (t rs -t rn )/(t rn + 1 -t rn ) + 100 × n. Literature RIs were obtained from http://webbook.nist.gov.

DNA extraction
DNA was isolated from 16 populations (Table 1) using a slightly modified procedure previously described as protocol D (Dodos et al. 2014), with the following modifications: the use of polyvinylpolypyrrolidone (PVPP, AppliChem Gmbh, Germany) instead of PVP 10 in the extraction buffer, and the use of 450 μL of 5 M NaCl in the isopropanol step. DNA purity and quantity was determined using a LambdaBio (PerkinElmer, UK) spectrophotometer. Samples were stored at -80 °C until further analysis.
Macrogen Europe B.V. performed fragment analysis using ABI 3730XLs capillary DNA analyser. The obtained electropherograms were analysed using Peak Scanner Software 2 (Applied Biosystems). Only fragments between 50 and 600 bps with a relative height above 100 were scored. Polymorphic bands were scored as present (1) or absent (0) and formed a binary matrix. A present/absent matrix of fragments obtained by all three selective primers was exported and processed in the R statistics software package (R version 3.5.1), where unreliable characters were removed and mistyping error was estimated using RawGeno 2.0 (Arrigo et al. 2012). The processed matrix was then used in statistical analysis.

Environmental data
Bioclimatic data based on the 30-year average temperature, precipitation, solar radiation and wind speed of all 18 studied localities were taken from the WorldClim 2.0 set of global climate layers with the precision of ~1 sq km (Fick and hijmans 2017). Data were analysed using QGIS software (QGIS Development Team 2015). Geological substratum information was obtained from the International Geological Map of Europe and Adjacent Areas (Asch 2005). The soil type and ph were retrieved from ISRIC (hengl et al. 2014).

Statistical analysis
Standard statistics (mean, standard deviation and distribution) were used to study data before performing univariate  (2) 17,173 l a mixed sample of ten individuals; sub. -substrate type; l -limestone; D -dolomite; P -phyllite; c -conglomerate; g -gneiss; u -undifferentiated; incl -inclination (in %); alt -altitude, m.a.s.l.; eo -number of individuals for essential oil analysis; aflP -number of individuals for aflP analysis, numbers in parentheses represent the number of individuals used in statistical analysis; Beou -number of vouchers in Beou herbarium; longitude and latitude in degrees (Analysis of variance, ANOVA) and multivariate analyses (principal components analysis, PCA; discriminant analysis, DA; hierarchical cluster analysis, hCA; multivariate analysis of variance, MANOVA). For all statistical analyses, components present above 0.1% on average that were not correlated (-0.7 < R > 0.7) and had normal distribution were used. Simple linear correlation and multivariate correlation (Mantel and partial Mantel tests) were used to analyse the correlation of EO components with environmental data. In the Mantel test, squared Mahalanobis' distances were used for EO data, Euclidean distances for log-transformed bioclimatic data and geographical distances based on GPS coordinates were used. All statistical analyses were performed using PAST 4.07 (hammer et al. 2001).
For AFLP data, AMOVA (Analysis of Molecular Variance) was calculated using GenoDive v.3.0 (Meirmans 2020), with 9999 permutations. Population structure based on a squared Euclidean distance matrix using Fst-analogue independent of ploidy level and breeding system (Ronfort et al. 1998;Meirmans and Liu 2018) was analysed. Bayesian model-based estimation of population structure was performed using Structure version 2.3.4 (Pritchard et al. 2000). The analyses were performed under the admixture model assuming independent allele frequencies and using a burn-in period of 10,000 followed by 5000 Markov Chain Monte Carlo, and the most likely value of K was determined according to Evanno et al. (2005). Principal coordinate analysis (PCoA) was performed to illustrate overall similarity among individuals using GenAlEx 6.5 (Peakall and Smouse 2006). PCoA was inferred from Nei's pairwise genetic distance (Nei 1978) between all pairs of AFLP phenotypes. The Mantel test was carried out to assess the correlation between genetic (Nei's genetic distances of populations) and geographic distances.

Environmental factors
Previous studies of EO in junipers have shown their chemophenetic significance (Adams 2000;Rajčević et al. 2013Rajčević et al. , 2015Rajčević et al. , 2018Rajčević et al. , 2020. however, the initial analysis of the literature data showed variation in EO composition based on the locality (see Section "EO composition"). To exclude the potential influence of environmental factors on EO composition in the studied populations, 102 environmental variables were considered. Analyses were carried out using raw data (376 individual samples) and population data (mean population values, including Mt. Orjen).
Based on the bioclimatic parameters, localities can be separated into three distinct groups (Supplemental Figure S1). Localities from Slovenia and northern Croatia share similar bioclimatic characteristics -intermediate temperatures and more precipitation than the other two groups. The second and third groups have lower precipitation in common but differ in temperatures. The populations from Serbia, Romania and Bulgaria grow in the coldest region of the three. In contrast, all other populations (Greece, Montenegro, Bosnia and herzegovina, and Northern Macedonia) live in a warmer climate. In addition, the localities differed in the dry periods. In addition to bioclimatic parameters, geographic position, inclination, exposition, wind speed, solar radiation, geological substrate, soil type and soil ph were also analysed. Populations from Mt. Stara, Mt. Velebit and Mt. Biokovo were characterised with the highest wind speeds (more than 5 m/s on average), Mt. Risnjak, Mt. Snežnik and Mt. Orjen intermediate (4-5 m/s on average), while all others had lower wind speeds (2-4 m/s) throughout the year (Supplemental Figure S2a). On the other hand, solar radiation was more-or-less similar throughout the study area, with a gradual decline correlated with higher latitudes (Supplemental Figure S2b).
The majority of the populations had limestone as substrate, with several growing on phyllite, dolomite and gneiss (Table  1). In almost all cases, the soil type was cambisol or a cambisol/podzol mixture. however, there was no separation of populations related to these parameters, e.g. populations from Mt. Kopaonik, Mt. Stara, Mt. Suva grow on different mother rocks but could not be differentiated based on the EO composition. At the same time, populations growing on the same soil type showed a different distribution of chemotypes.

AFLP analysis
From the 16 populations, 90 individuals were chosen for AFLP analysis. To test the reproducibility of the AFLP markers, we repeated the entire AFLP procedure twice with 23 randomly chosen individuals. Forty-three individuals with the strongest peaks in the electropherograms were selected for statistical analysis. The mistyping error was low -0.01%, as calculated by the RawGeno 2.0. The three selective AFLP primer pairs generated a total of 642 fragments. Of these, 641 (99.8%) were polymorphic. Both polymorphic and monomorphic fragments were used for the calculations. The number of monomorphic loci in each studied population varied from 8.3% in Mt. Galičica to 38.6% and 38.9% in Mt. Ţarcu and Mt. Stara, respectively. Most of the loci (86.6%) were present in all regions, while only five were private for regions, and 80 were present in two out of three regions. 14.5% of loci were consistently present in all analysed populations, and only 4.3% were present in 95% of the individuals.
AMOVA showed that most of the genetic variability was within the populations (91.7%), whereas the rest was between populations and the three bioclimatic regions (Table 2). These results are in concordance with the previous studies carried on populations from Britain, northwestern Europe and central Europe (Merwe et al. 2000;Michalczyk et al. 2010;Vanden-Broeck et al. 2011), where authors found high intrapopulation diversity and low differentiation between populations.
The results of AMOVA were further confirmed using the PCoA with Nei's distances. The first two eigenvectors explained 23.1% of the total variability. The PCoA scatter plot showed the formation of three clusters (Figure 1), though these clusters did not correlate with geographic regions, which was also confirmed by the Mantel test (R = -0.12, p = 0.79). Even though there is no separation of populations, these results show that most of the analysed individuals from some populations tend to cluster closer together, regardless of their sex or EO composition. Structure analysis also confirms this. K-cluster analysis plateaued quickly, and based on the Bayesian information criterion, two K-values had the highest probability (Figure 2). In both cases, most individuals have mixed ancestry, suggesting recent or constant gene flow between these populations.

EO composition
In the EO obtained from 376 individuals of J. communis var. saxatilis, 156 components were detected and identified, representing on average 99.0 ± 0.9% of the total EO composition    Table S1). The number of compounds per sample ranged from 51 to 101. The chemical composition was dominated by monoterpenes, particularly monoterpene hydrocarbons: up to 72.9% on average in the Mt. Kopaonik population. Only in three populations (Mt. Risnjak, Mt. Snežnik and the Julian Alps) were sesquiterpenes more abundant, accounting for 47.1 to 59.7% of the total oil composition. Diterpenes were present in low amounts (≤6%) in most populations, except in the two easternmost populations, Mt. Balkan and Mt. Țarcu, with 9.5% and 18.2%, respectively.
The majority of the populations had high quantities of α-pinene, ranging from 16.2 to 49.0%. On the other hand, sabinene was present in low abundance (<5%) in eight populations and in medium-to-high amounts (>5%) in ten. The population from Mt. Kopaonik had the highest amount of sabinene, 22.0%. δ-3-Carene was present in five populations in medium-to-high abundances (the highest was 14.7% in Mt. Risnjak). Germacrene D was present in low-to-medium abundances in all populations (2.6%-12.0%, highest in Mt. Snežnik).
Analysis of EOs has shown that individuals that grow only a few meters apart can have different dominant compounds, thus corroborating the chemophenetic significance of EOs. The previous investigation of J. communis var. communis in the north-western Balkans (Rajčević et al. 2020) showed significantly different EO profiles -these EOs were strongly dominated by monoterpenes (54%-75%), with a strong dominance of monoterpene hydrocarbons (48%-58%), in contrast with the populations of var. saxatilis from the same region (monoterpenes 37%-55%, of which almost all were monoterpene hydrocarbons), so we compared the EOs from the present research only with EO compostion of var. saxatilis from the literature. Previous investigations of Juniperus communis var. saxatilis have revealed several potential chemotypes throughout its range, with α-pinene (Italy, Greece, Bulgaria, Mongolia and Korea), sabinene (Portugal, Switzerland, Corsica, Serbia and Iran) and limonene (Italy, Corsica and Sardinia) being the most common (Table 3). δ3-Carene had a high abundance only in Greece, though it was not the most abundant compound in the EO. however, most of these studies dealt with plant material collected from several individuals and processed in bulk, thus representing the average EO for those collected individuals or a single individual (Bulgaria, Serbia). Taking into account the detected variability, one must be very careful when drawing conclusions. Only Ottavioli et al. (2009) studied intraindividual variability and identified two chemotypes: sabinene and limonene. Interestingly, the highest limonene content in the present work was found in a single individual from Mt. Durmitor (14.9%), while all the others had this compound at fairly low concentrations (around 2%), clearly distinguishing the populations in Corsica from those in the Balkans.

EO variability
Fifty-five compounds that were present above 0.1% on average and were not correlated (-0.7 > R < 0.7) were used for further statistical analyses. The first two eigenvectors of PCA explained 81.9% of the total variability, indicating separation of individuals based on three compounds: α-pinene, sabinene and δ-3-carene (not shown). however, there was no apparent separation of populations, geographic regions or sex in this analysis, suggesting high intraspecific and intrapopulation variability.
MANOVA and DA were also performed to assess the differentiation of populations. In MANOVA, the mixed sample from Mt. Orjen was omitted, while the population from Mt. Snežnik was merged with the population from Mt. Risnjak due to their extreme proximity and a small number of samples. The MANOVA showed a statistically significant differentiation of the populations (F = 11.3, p = 8.9 ×·10 −8 ) based on 55 components, though the post hoc test (pairwise hotelling's tests) showed the formation of some groups, i.e. the MANOVA was not able to differentiate populations from Mt. Kopaonik, Mt. Suva and Mt. Stara from each other, nor the populations from Biokovo and Mt. Prenj. Their geographic proximity can easily explain this. however, MANOVA was also unable to differentiate three populations: Mt. Galičica, Mt. Snežnik + Risnjak and Mt. Olympus, even though they are geographically separated. DA showed the formation of three distinct clusters based on the presence of α-pinene, sabinene and δ-3-carene (Figure 3a). The populations from Mt. Prenj, Mt. Durmitor, Biokovo, Mt. Velebit and Mt. Bistra were grouped based on the higher abundance of α-pinene. In the post hoc test, Mt. Orjen was placed within this group. The second group included populations with a higher amount of sabinene: Mt. Komovi (Montenegro) and all populations from Serbia, Romania and Bulgaria. The third cluster is separated from all others based on the highest abundance of δ-3-carene (and the lowest of α-pinene and sabinene): Mt. Snežnik, Julian Alps, Mt. Galičica, Mt. Olympus and Mt. Smolikas. These populations also had higher abundances of different germacrene compounds (Supplemental Table S1). Based on the EO composition, 91.5% of individuals were correctly assigned to their predefined groups (i.e. populations) in the confusion matrix, suggesting some separation between them.
Since PCA and DA strongly suggested the presence of, at least, three chemotypes, individuals were assigned to one of three groups based on the relative main components from PCA and DA and the abundance of α-pinene, sabinene and δ-3-carene -putative chemotypes. ANOVA, MANOVA and DA were performed to assess the distinction of these chemotypes. ANOVA showed a statistically significant separation of 40 of the 55 compounds (Supplemental Table S2) and all of the chemical classes, while MANOVA showed a statistically significant separation of all three chemotypes (p = 3.9 × 10 −165 ), including post hoc tests (pairwise hotelling's tests). DA (Figure 3b) showed almost complete separation of the three putative chemotypes, with 98.4% of individuals correctly classified to their assigned group (chemotype). Only six individuals were re-classified, five of them from the δ-3-carene group. It is worth mentioning, that the δ-3-carene group had a quite variable composition that was not always characterised with the dominance of δ-3-carene. This designation was assigned due to the most significant component according to PCA and DA analyses. This chemotype was characterised by almost equal amounts of sesquiterpenes and monoterpenes, namely due to high abundances of different germacrenes (e.g. germacrene D, germacrene B, germacrene A and germacrene D-4-ol), so it could also be designated as the germacrene chemotype.
In their work on the correlation between EO profile and grazing by herbivores, Markó et al. (2008) also found these three chemotypes in J. communis var. communis. In this research, authors found that the individuals that produced a high amount of δ-3-carene were undamaged by the local herbivores. So the presence of a greater number of individuals with this chemotype in some areas suggests that these populations might have been used for sheep or another herbivore grazing in the past since no such activity was detected in the analysed populations, or was noted by the authors during the collection of the plant material.

Correlation of EO with environmental variables
Using raw data, no correlation was found between EO composition and site bioclimatic characteristics, which was expected given the high variability of all EO components among populations (data not shown). however, when using the mean population values for EO components, the Mantel test found a significant but low correlation (r = 0.41, p < 0.01) with the bioclimatic parameters of the localities, suggesting some influence of climate on the distribution of chemotypes in populations across the Balkans. Linear correlation test (Supplemental Table S3) showed a significantly high correlation of some of the compounds with the bioclimatic data, e.g. α-pinene was strongly positively correlated with mean annual temperature, whereas sabinene and myrcene showed Figure 3. Da scatter plot (a) populations as groups, with three main vectors drawn, 1 -the Julian alps, 2 -mt. snežnik, 3 -mt. risnjak,4 -mt. Velebit, 5 -Biokovo,6 -mt. Durmitor, 7 -mt. Komovi, 8 -mt. orjen, 9 -mt. Prenj, 10 -mt. Kopaonik, 11 -mt. suva, 12 -mt. stara, 13 -mt. Ţarcu, 14 -mt. Balkan, 15 -mt. Bistra, 16 -mt. galičica, 17 -mt. olympus, 18 -mt. smolikas. for population details cf. high negative correlations with precipitation. Groups of compounds in EO also showed a correlation with some of the bioclimatic parameters. The abundance of monoterpene hydrocarbons showed a high positive correlation with the temperature variability and precipitation seasonality, i.e. the greater the temperature and precipitation variability during the year, the greater the abundance of monoterpene hydrocarbons in EO composition. On the other hand, precipitation of the driest month and quarter showed a high negative correlation with the abundance of monoterpene hydrocarbons. Of the remaining terpene groups, only diterpene hydrocarbons showed a significant and high correlation with precipitation seasonality. All other groups were not correlated with the bioclimatic parameters.
Mantel test showed a moderate correlation between geographic region and EO composition (R = 0.47, p < 0.01), while the linear correlation (Supplemental Table S3) was high and positive for sabinene and longitude, as well as for several sesquiterpene compounds and longitude. Although only monoterpene hydrocarbons showed a high positive correlation with latitude, all compound groups showed a moderate-to-high correlation with longitude -monoterpenes and diterpenes positive, sesquiterpenes a negative correlation, i.e. the further east the population grows, the lesser amount of sesquiterpenes are produced on average. In the linear correlation test, inclination, exposition and altitude did not show any significant correlation with any of the components of EO, with the exception of altitude with sabinene and trans-totarol (R = 0.53 and 0.60, respectively, p < 0.01). Several compounds had a high correlation (-0.6 < R > 0.6, p < 0.01) with the soil ph: β-pinene (positive), germacrene B, germacrene D-4-ol, and abieta-8,12-diene (negative).
Even though environmental factors did not show correlation with the EO composition at the individual level, the distribution of the chemotypes is correlated with the environmental factors, of which precipitation and temperature are the most important ones (Figure 4). Based on the bioclimatic parameters, two main regions can be distinguished -a region with higher precipitation in the northwestern Balkans and a drier climate in the rest of this region. Furthermore, the drier areas of the Balkans can be divided into warmer and colder regions based on temperatures. In the north-western Balkan region, δ3-carene chemotype was most prevalent, with several individuals having the sabinene chemotype. Interestingly, populations of typical variety from that region (var. communis), often growing at the same locality with J. communis var. saxatilis, showed domination of sabinene chemotype with relatively high abundance of limonene and extremely low abundance of δ-3-carene (Rajčević et al. 2020).
In addition, sesquiterpene content was the highest in these populations, which was also confirmed in correlation tests, which was not the case with var. communis. however, the population from Mt. Velebit has a chemotype distribution more similar to that of other coastal mountains in the Adriatic. This could be explained by the high wind speeds in the locality, which contribute to higher local aridity, regardless of precipitation. Also, this locality had the highest precipitation in the colder part of the year, so most of the precipitation consisted of snow. In contrast, α-pinene and sabinene chemotypes were more common in the drier areas in the Balkans. No sabinene chemotypes were detected in two coastal mountains (Biokovo and Orjen) and two adjacent mountains (Mt. Prenj and Mt. Durmitor). On the other hand, this chemotype (sabinene) was more common in the central and eastern Balkans. Mt. Komovi, though relatively close to the coast, is surrounded by a massive mountain range, contributing to a much drier climate.

Conclusion
EO analysis of J. communis var. saxatilis obtained from 376 individuals from 18 populations shows high interpopulation variability. This variability is not correlated with sex or the studied environmental factors but is in concordance with the Figure 4. spatial distribution of analysed populations in south-eastern europe. Pie charts show the ratio of detected chemotypes for each locality. orange -δ3-carene, red -sabinene, blue -α-pinene. for population details, cf. table 4. molecular data. J. communis var. saxatilis has a fragmented distribution in the Balkans, occurring only at high altitudes. Nevertheless, the molecular data showed high variability within populations and low differentiation between populations. This is in agreement with most of the studies dealing with long-lived woody species, which tend to show abundant gene flow and long generation times. Wind pollination and seed dispersal by birds of this species also contribute to high gene flow between geographically separated populations (Jimenez et al. 2017).
Even though there was no correlation between EO composition and the environmental variables studied over extended periods, it appears that natural selection made certain chemotypes more common despite the high gene flow and genetic diversity of this species. The high genetic variability in the region and low population differentiation indicate that there is no immediate threat to populations of prostrate junipers. however, as the climate in the region becomes drier, one might expect a shift in the abundance of the different chemotypes towards chemotypes containing more sabinene and a decline in δ3-carene chemotypes. This will certainly affect variability and endanger better survival under future environmental changes. In addition, all three compounds give different aromas to the extracts, α-pinene woody, sabinene spicy and δ3-carene sweet (Opdyke 1979). They also have different medicinal properties. All this can potentially lead to problems for the local farmers/collectors relying on certain chemotypes for the spirit industry and, also, for the use of these plants in traditional medicine, which is especially important in rural mountainous areas. To assess the potential adverse effects of this change additional studies are needed.