Patterns of Genetic Variation of Nepeta nuda L. from the Central Balkans: Understanding Drivers of Chemical Diversity

Nepeta nuda L., a notable medicinal species in the tradition of the Balkan region, is a rich source of bioactive iridoids and phenolics previously described as high-resolution taxonomical classifiers for the genus Nepeta. However, their potential in investigating intra-species differentiation is here described for the first time. The aim was to recognize the sources of natural chemical diversity and their association with the genetic variability both within and among N. nuda populations in the Central Balkans. Chemical diversity was assessed from methanol extracts and essential oils through untargeted and targeted metabolomics using state-of-the-art analytical tools, covering a broad spectrum of compounds that represent the N. nuda metabolome. We found that chemodiversity primarily resides within populations of N. nuda, and similar results were obtained at the DNA level using microsatellite markers. The low genetic and chemical differentiation of the studied N. nuda populations implies that their metabolomic profiles may be less influenced by geographic distance and variable environmental conditions within the Central Balkans, as they are under the pivotal control of their genetic backgrounds. Screening the distribution of the major bioactive compounds belonging to phenolics (phenolic acids and flavonoids) and iridoids (both aglycones and glycosylated forms), within and among N. nuda populations, is able to guarantee mass spectrometry-based tools for the selection of elite representative genotypes with practical importance. The knowledge acquired will allow us to delve deeper into the molecular background of N. nuda chemical diversity, which is the course of our further work.


Introduction
Nepeta nuda L. (fam.Lamiaceae) is a widely distributed Eurasian species, traditionally used in folk medicine of the Balkan nations [1], along with other congeneric taxa.This remarkable plant can be found in forest clearings and meadows, at montane and subalpine altitudes up to 2100 m [2], and often in ruderal habitats alongside roads where its robust habitus, represented by colorful long-lasting flower spikes, dominates other species in the surrounding (personal observation).N. nuda is traditionally used as a remedy for various ailments such as gastrointestinal problems and respiratory diseases, as well as for treating wounds [3].Numerous studies have shown that extracts of N. nuda have antiviral [4], antioxidant [5], and antibacterial [6,7] properties.In addition, essential oils of N. nuda have shown phytotoxic and herbicidal effects [4,8], which make this species potentially useful as a novel source of potent biopesticides.The various biological activities of its extracts and essential oils could be attributed to the two dominant classes of specialized metabolites: terpenoids (mainly monoterpenes and sesquiterpenes) and phenolics (primarily flavonoids and phenolic acids).The most abundant monoterpenes among Nepeta species occur either as aglycones or glycosides of iridoids.A characteristic representative of the group of aglycones is nepetalactone, a bicyclic monoterpenoid exclusively found in Nepeta species that occurs in the form of different diastereoisomers but also as hydrogenated derivatives such as 5,9-dehydronepetalactone or nepetalic acid.Iridoid glycosides, featured by the presence of a sugar component attached to an iridoid scaffold, are 1,5,9-epi-deoxyloganic acid, nepetaside, and nepetariaside (literature survey is summarized in Supplementary Table S1).Chemical composition of N. nuda can be influenced not only by the geographical origin and harvesting time but also by environmental parameters such as temperature, precipitation, and solar radiation [3,[9][10][11].Moreover, variable qualitative and quantitative composition of nepetalactone stereoisomers can be found in N. nuda essential oils (EOs) of different origins (Supplementary Table S1), which indicates that the metabolome of this species might be highly susceptible to the influence of both endogenous and environmental factors.According to Aćimović et al. [3], nepetalactones represent minor components of the essential oils of N. nuda collected in Serbia, while in populations from Greece, Turkey, and Iran, various nepetalactone stereoisomers constitute a large proportion of their essential oils [4,12,13].The most characteristic phenolic constituents of N. nuda are rosmarinic, caffeic, chlorogenic, and ferulic acids, as well as aesculin [10,14,15].However, chemical profiling of non-volatile compounds such as phenolics or iridoid glycosides has not been extensively performed, as most analyses were focused on the characterization of N. nuda essential oils composition (Supplementary Table S1).
N. nuda has a wide geographical range [16] and inhabits regions with different climate conditions, soil types, and altitudes, while various biotic and abiotic factors affect its fitness and biomass productivity.Increasing evidence shows that maintaining genetic diversity within natural populations of various plant species can maximize their potential to withstand and adapt to environmental fluctuations [17], although it is difficult to predict the novel selection pressures to which populations will be exposed.The evaluation of the genetic diversity of natural populations of N. nuda is directed into providing the evidence on their differentiation magnitude and predicting the likelihood of their survival within the contemporary environmental pressures such as global climate change.Addressing a great challenge of mankind to increase agricultural productivity while preserving and increasing biodiversity [18], we bring into focus the relationships between the genetic and chemical diversity patterns of N. nuda populations from the Central Balkans, thus setting the stage for the identification and selection of lines featuring key traits for maximizing the biomass and secondary metabolite production of this medicinal plant.Understanding the reciprocal relationships between the genotype and phenotypic traits in this species will allow us to extend our findings to other congeneric species.Recognition of highly and less-productive populations and genotypes of N. nuda will facilitate our quest to explain the molecular background of the chemical diversity of iridoids within the genus Nepeta, which is the course of our future work.To the best of our knowledge, no studies dealing with N. nuda intra-species chemodiversity with reference to its genetic background in wild populations have previously been conducted, and our intention was to fill these knowledge gaps by adopting the combination of EST-SSR and SSR molecular markers in parallel with state-of-the-art metabolomics tools which are able to encompass different components of the N. nuda metabolome.Codominant markers have not been previously employed to estimate the genetic variations of any Nepeta sp. at the population level through the estimation of heterozygosity, allelic fixation index, variance fractioning, or population structuring.

Results and Discussion
Qualitative and quantitative profiling of specialized metabolites, while capturing genotype-dependent variations, enabled us to comprehensively acquire, analyze, and adequately interpret the overall chemical diversity of N. nuda in the Central Balkans.State-of-the-art analytical techniques (LC-MS and GC-MS) that cover different parts of the metabolomes of the targeted taxon are well-certified tools to comprehensively assess the underlying chemical diversity and complexity and to aid in predicting differences at interand intra-population level [19].As chemodiversity has its foundations in genetic diversity, we utilized a set of EST-SSR and SSR markers to produce informative DNA data for the reconstruction of the genetic background of N. nuda populations which, to our knowledge, has not been undertaken before.

UHPLC-ESI-QToF-MS Identification of N. nuda Leaf Metabolites
The subject of this research was the detailed non-targeted characterization of metabolites accumulated in the leaves of N. nuda, inhabiting eleven locations in the Central Balkans (Figure 1).The focus was on phenolic and iridoid compounds, which are the two most significant groups of bioactive metabolites in the Nepeta species [19,20].Data on previous records about iridoids and phenolics previously reported for N. nuda are summarized in Supplementary Table S1.Using an UHPLC-ESI-QToF-MS instrument operating in both positive and negative modes, seventy-eight compounds were identified, which belong to hydroxybenzoic (ten compounds) and hydroxycinnamic (nineteen compounds) acid derivatives, iridoid glycosides (nine compounds), iridoid aglycones (sixteen compounds), flavonoid glycosides (two compounds), and flavonoid aglycones (six compounds), as well as sixteen compounds from other groups of metabolites.

Results and Discussion
Qualitative and quantitative profiling of specialized metabolites, while capturing genotype-dependent variations, enabled us to comprehensively acquire, analyze, and adequately interpret the overall chemical diversity of N. nuda in the Central Balkans.Stateof-the-art analytical techniques (LC-MS and GC-MS) that cover different parts of the metabolomes of the targeted taxon are well-certified tools to comprehensively assess the underlying chemical diversity and complexity and to aid in predicting differences at interand intra-population level [19].As chemodiversity has its foundations in genetic diversity, we utilized a set of EST-SSR and SSR markers to produce informative DNA data for the reconstruction of the genetic background of N. nuda populations which, to our knowledge, has not been undertaken before.

UHPLC-ESI-QToF-MS Identification of N. nuda Leaf Metabolites
The subject of this research was the detailed non-targeted characterization of metabolites accumulated in the leaves of N. nuda, inhabiting eleven locations in the Central Balkans (Figure 1).The focus was on phenolic and iridoid compounds, which are the two most significant groups of bioactive metabolites in the Nepeta species [19,20].Data on previous records about iridoids and phenolics previously reported for N. nuda are summarized in Supplementary Table S1.Using an UHPLC-ESI-QToF-MS instrument operating in both positive and negative modes, seventy-eight compounds were identified, which belong to hydroxybenzoic (ten compounds) and hydroxycinnamic (nineteen compounds) acid derivatives, iridoid glycosides (nine compounds), iridoid aglycones (sixteen compounds), flavonoid glycosides (two compounds), and flavonoid aglycones (six compounds), as well as sixteen compounds from other groups of metabolites.UHPLC-ESI-QToF-MS data on metabolites identified in the leaf extracts of eleven N. nuda accessions are summarized in Table 1.The base peak chromatograms, in negative and positive ionization modes, of 11 different Nepeta nuda leaf extracts are depicted in Supplementary Figure S1.References related to the previous identification of these compounds in different Nepeta species are listed as well.By searching the available literature, we found that seven compounds (Supplementary Figures S2-S5) listed in Table 1 were identified for the first time in any Nepeta species.Supplementary Table S2 summarizes the peak areas of the identified compounds in the eleven investigated extracts.
From the group of ten hydroxybenzoic acid derivatives, four compounds were identified as free acids (2, 3, 6, and 9), while the other six, especially abundant, were found to be hexosyl derivatives, giving a specific fragment ion resulting from the loss of 162 Da (hexosyl residue).The fragmentation pathway of compound 5 (vanillic acid hexoside), for which we found no data, confirming that it has not previously been identified in any Nepeta species, is shown in Supplementary Figure S2A.
All nineteen compounds from the group of hydroxycinnamic acids, except for coumarin umbelliferone (compound 16), are actually different derivatives of caffeic acid.The largest number of these derivatives were identified as esters with tartaric, tartronic, glycolic, or malic acids, which are otherwise common for Nepeta species [24][25][26]29].Compound 25, which was not previously identified in Nepeta species, with the molecular ion at 309 m/z, showed specific MS 2 ions corresponding to deprotonated ferulic acid (193 m/z) and its fragments (178, 147, 134, and 117 m/z) (Supplementary Figure S2B).Nine derivatives of iridoid glycosides were identified in negative ionization mode (Figure 2A).Some compounds showed [M-H] − as a molecular ion, while in three iridoid glycosides (32, 34, and 35), an adduct with formic acid ([M +HCOOH-H] − ) was observed, giving a shift in the molecular ion mass of 46 Da.The most abundant compound, compared by peak area, was nepetanudoside A (compound 35).This compound, typical for N. nuda [33], was detected in all but one population (Straža) in our research.One of the derivatives of iridoid glycoside (compound 33, 5-deoxylamiol) has not been previously reported in any Nepeta species but is typical for Lamium species from the same family [52].The proposed detailed fragmentation pathway of this compound is shown in Supplementary Figure S3.
Iridoid aglycones were mostly detected in the positive ionization mode due to their polarity, except for 7-deoxyloganetin and nepetalic acid (compounds 49 and 50), which gave more abundant peaks in the negative ionization mode.Compound 49 was identified in only one population (Brodica), while compound 50 was detected in all populations.Out of a total of sixteen identified compounds, the presence of fourteen was confirmed by comparison with the available literature on Nepeta species (Table 1), while genipin and loganetin (compounds 39 and 42) were identified through the evaluation of their MS data.Supplementary Figure S4A represents the detailed fragmentation pattern in the positive ionization mode of genipin (compound 39), which was detected in three populations: Selište, Donji Krivodol, and Gornje Selo.Compound 39 showed MS 2 the base peak at 103 m/z in the positive ionization mode, which was in accordance with literature data [53].The iridoid loganetin (compound 42), which was not previously reported as common for the genus Nepeta, was detected in almost all populations (except for Straža) in considerable amounts.Its chemical structure and explanation of the formation of all MS 2 fragments is shown in Supplementary Figure S4B.
The vast majority of Nepeta species studied by Jamzad et al. [43] contains flavonetype flavonoids on the leaf surface.In the case of flavonoid glycosides, only two luteolin derivatives (compounds 55 and 56) were identified in our study, and both were previously reported in catnip [24,27].Our results also show the presence of six flavonoid aglycones belonging to this subgroup of flavonoids, while derivatives from other subgroups were not identified.By comparison of the peak areas of these compounds, it can be noticed that salvigenin (compound 62) was very abundant flavone.This compound was previously reported to be present only in the aerial parts of N. transcaucasica [44], in N. nuda, N. assurgens, N. asterotricha, N. crispa, N. pungens, N. saturejoides [43], and N. baytopii [29].
Ungrouped metabolites do not belong to any of the aforementioned classes.Two fatty acids (compounds 69 and 73) and four carboxylic acids (compounds 70, 72, 74, and 78) specific to the genus Nepeta were recorded (Table 1).Compound 75, with the molecular ion at 201 m/z, was identified as argolic acid A. It gave the MS 2 base peak at 139 m/z, resulting from the loss of H 2 O and CO 2 (62 Da).This cyclopentane-monoterpene diol was named after the species from which it was first isolated, N. argolica subsp.argolica [51].In our study, two new derivatives of argolic acid A were found and marked as argolic acid A rhamnoside and argolic acid A methyl ether rhamnoside (compounds 71 and 76, respectively).Their structures and fragmentation pathways are proposed in Supplementary Figure S5.To visualize relations in metabolic fingerprints of eleven N. nuda populations, HCA was performed based on the Spearman's method of cluster agglomeration.Population Donji Krivodol is clustered separately from other populations on the HCA plot (Figure 2B), owing to the accumulation of high amounts of caffeic acid and its derivatives.The remaining populations are separated into two major subclusters, the first one containing only the Straža population, characterized by high amounts of nepetalactone-like compounds.Populations belonging to the second subcluster are separated into two groups, the first one containing the populations Vlasina and Židilje and the second containing the remaining populations.PCA analysis segregated the populations Janjska reka, Gornje Selo, Brodica, and Donji Krivodol from the remaining ones along the PC1, which explains 29.78% of the total variability (Figure 2C).The major contributors to the diversification along PC1 are nepetanudoside A, caffeoyltartonic acid, 7-deoxyloganetic acid, salvigenin, and 1,5,9-epideoxyloganic acid.Nepetanudoside A, caffeoyltartronic acid, and caffeoylglycolic acid are the major factors behind the diversification along PC2, which explains 19.26% of the total variability and mostly affects the population Vlasina.PC1 and PC2 cumulatively explain 49.04% of the variability.

GC-MS Non-Targeted Metabolomics of N. nuda EOs
In total, 30 terpene compounds were identified in the EOs of N. nuda, adopting the GC-MS analysis (Table 2).
These compounds belong to the groups of monoterpenes (6 monoterpenes and 6 oxygenated monoterpenoids) and sesquiterpenes (16 sesquiterpenes and 2 oxygenated sesquiterpenoids).The representative GC-MS TIC chromatograms are shown in Figure 3A.Among monoterpenoids, 1,8-cineole (7), α-pinene (1), and β-pinene (3) predominated.From the group of monoterpenoid iridoids, two nepetalactone stereoisomers were detected: trans,trans- (11) and cis,trans-nepetalactone (12), albeit in low amounts.Interestingly, 1,5,9 -dehydronepetalactone and dihydronepetalactones, as well as their derivatives, were not recorded in the EOs of N. nuda through the GC-MS analysis.The same applies for the groups of diterpenes and triterpenes.This can be ascribed, at least partially, to the hydrodistillation method for the EO isolation but also to the HS sampling procedure, which includes the incubation of samples at 40 • C prior to sampling and injection.This temperature might be too low to induce the evaporation of these compounds, which might generally result in a lower abundance of less volatile compounds.beta-Caryophyllene (17) and Germacrene D (23) were detected as the major sesquiterpenoids in the majority of analyzed N. nuda EOs, while α-humulene (22) and β-bisabolene (24) were also present in significant amounts.All the analyzed EOs of eleven N. nuda populations displayed a similar chemotype characterized by the predominance of 1,8-cineole, germacrene D, and beta-caryophyllene.As the major constituents of N. nuda EOs, other studies reported caryophyllene oxide [54], camphor and 1,8-cineole [55], caryophyllene [56], nepetalactone [4,57,58], 1,8-cineole and 4aα,7β,7aα-nepetalactone [59], nepetalactone and germacrene [6], 1,8-cineole [60], and 1,8cineole and a mixture of nepetalactones and germacrene-D [61].Previously reported data on the composition of terpenes in N. nuda EOs are shown in Supplementary Table S1.These differences might arise from many factors, including plant origin and growth conditions, growth stage and plant parts analyzed, the subspecies or a cultivar in question, EO isolation procedure, and GC-MS conditions.When the composition of EOs of N. nuda sampled at Mt Rtanj (Serbia) was tracked through three successive years, the content and ratio of EO components were found to be slightly influenced by weather conditions such as temperature, precipitation, and insolation [62].The same authors proposed four chemotypes of N. nuda EOs: (1) nepetalactone; (2) 1,8-cineole; (3) mixed (nepetalactone+1,8-cineole+germacrene D); and (4) nonspecific chemotypes.Based on this principle, the Rtanj population within the present study, rich in 1,8-cineole, belongs to the chemotype 2. A previous study [11] revealed 4aα,7β,7aα-nepetalactone, 1,8-cineole, and germacrene D as the most abundant volatile compounds in N. nuda plants grown under either in vitro or ex vitro conditions.Caryophyllene, β-ocimene, bicyclogermacrene, β-pinene, myrcene, and humulene were also found to be abundant.Although a previous study reports the comparative analysis of EO composition in four N. nuda populations from Iran [13], the present study represents the first successful characterization of interpopulation variability in the qualitative EO composition of this species.

Targeted qqqMS Profiling of Major Iridoids and Phenolics in N. nuda Methanol Extracts
To estimate inter-and intra-population phenotypic diversity in the quantitative composition of major phenolics and iridoids in N. nuda, and thus the productivity of individual genotypes, we further applied a targeted metabolomics approach.Selected phenolics (six compounds) and iridoids (three compounds) were previously adopted as high-resolution taxonomic classifiers to discriminate Nepeta species [20].The set selected for the present study comprised a total of 12 compounds from the group of iridoids (four compounds) and phenolics (eight compounds).Although some previous studies revealed that the metabolic profiles of phenolics and iridoid glycosides depend on variable growth conditions [10], this is the first report on intra-and inter-population comparative analysis of these metabolites in N. nuda.
Among the targeted phenolic compounds, rosmarinic acid was the most abundant in all analyzed populations (Figure 4A).This was in accordance with previous studies, which highlighted rosmarinic acid as the major phenolic compound in N. nuda [10,14,15,20,63] and other Nepeta species [19,20,64,65].Samples from Janjska reka displayed the highest amounts of this phenolic acid (6698 µg 100 mg −1 DW) followed by samples from Gornje Selo, Židilje, and Straža (Figure 4A).A similar trend was observed for caffeic acid, where Janjska reka displayed the highest amounts of this compound (~4 µg 100 mg −1 DW).Chlorogenic acid was the most abundant in samples from Selište, with average values reaching 0.13 µg 100 mg −1 DW.Among the analyzed flavonoids, astragalin was the most abundant, especially in samples from the population Židilje (~11.70 µg 100 mg −1 DW).Samples originating from Vinatovača were the richest in apigenin (0.21 µg 100 mg −1 DW) and luteolin (0.15 µg 100 mg −1 DW).As for rutin, the amount of this compound showed prominent variability within populations; thus, no significant differences in the average amounts were recorded between populations.
To evaluate the variability in the quantitative composition of iridoids and phenolics in N. nuda within the Central Balkans area, PCA was performed, adopting Ward's method of the data agglomeration (Figure 4B).The major contributors to the diversification of samples along the PC1, explaining 97.45% of the total variability, are cis,trans-nepetalactone, rosmarinic acid, 5,9-dehydronepetalactone, and 1,5,9-epi-deoxyloganic acid.As for the PC2, which contributes to the total variability with 2.07%, nepetanudoside A is the major contributor to the diversification along this component.According to the HCA plot (Figure 4C), constructed via Pearson's algorithm, population Vinatovača was segregated from the other populations, as it was the richest in 1,5,9-epi-deoxyloganic acid, 5,9-dehydronepetalactone, apigenin, and luteolin.The remaining populations were grouped into two subclusters, the first one contained populations Donji Krivodol, Gornje Selo, and Selište, while the second contained the remaining seven populations.Interestingly, the results display high intra-population variability, especially for Janjska reka, Gornje Selo, and Brodica.Although iridoids are probably the most potent bioactive compounds in Nepeta species, studies dealing with simultaneous analysis of iridoid aglycones and their glycosides are scarce, most probably due to different extractions procedures and analytical techniques required for their precise identification and quantification.Here, we adopted cutting-edge analytical instrumentation to analyze in parallel iridoid aglycones and glycosides in N. nuda accessions.When compared to other congeneric species, N. nuda can be described as a low producer of iridoid aglycones (nepetalactones and their derivatives) but a relatively significant producer of glycosylated iridoids.Interestingly, genotypes containing iridoid aglycones in trace amounts are also present among the analyzed populations, which might have important consequences on their performance in natural habitats in the way they chemically respond to biotic constraints.The observed diversity in iridoid content might be due to the genotype-dependent regulation of iridoid biosynthesis in N. nuda, which, in other Nepeta species, usually occurs at the transcriptional level and affects the expression of biosynthetic genes and the corresponding transcription factors [66,67].Our further work is aimed at deciphering the regulatory mechanisms that define the metabolic flux through the iridoid pathway branches leading to nepetalactone and its derivatives or to iridoid glycosides, among which nepetanudoside A and 1,5,9-epi-deoxyloganic acids predominate.This will be assisted by simultaneous metabolic profiling and gene co-expression analysis of the highly and less-productive N. nuda genotypes identified within the present study.The targeted metabolic profiling of major iridoids and phenolics indicated that some of individuals/genotypes were characterized by unique metabolic profiles, which may be selected for further profiling of expression of the genes related to the biosynthesis of iridoids and rosmarinic acid.This is the course of our further work.
On the other hand, N. nuda is a rich source of rosmarinic acid, previously reported to be the most abundant phenolic compound among Nepeta species.Within the present study, it was possible to identify N. nuda genotypes rich in rosmarinic acid, which will be used in our quest for the factors regulating the metabolism of this phenolic acid.Regulatory mechanisms most likely vary across different biosynthetic pathways, such as terpenoids and phenolics (Zhiponova, M. personal communication), and exploring the correlations between these two classes of metabolites, at the level of metabolomes, transcriptomes, and proteomes, may support the development of computational models that aim to link chemodiversity to genetic and evolutionary principles.

Parameters of Genetic Variation within and among the Populations
It has been well documented that genetic background plays a key role in defining the composition and overall yield of specialized metabolites in plants.However, genotype and environment interactions determine the overall composition of metabolites at multiple levels of biological organization.The aim of the current study was to develop a knowledge resource for N. nuda chemodiversity, with reference to the genetic background in wild populationstowards facilitating future utilization for commercial and pharmacological applications.
To elucidate the fundamentals of species' natural chemical diversity, we assessed the genetic structure of wild populations, adopting a set of EST-SSR molecular markers mined from a leaf transcriptome of N. nuda and a pair of previously published SSR markers.Microsatellites, being highly polymorphic codominant markers with a neutral evolutionary history, provide a basic knowledge of the species' genome variation through the estimation of the allelic fixation index, as well as the structure of populations, offering outstanding benefits because of their low cost and user-friendly nature [68].EST-SSRs markers, being transcriptome-derived, are considered less variable as compared to genomic microsatellites, but meticulous selection of loci to be used can result in a greater discriminating power in comparison to genomic SSRs [69].To the best of our knowledge, no prior attempts to characterize genetic diversity of any taxon belonging to the genus Nepeta employing microsatellites was made.Previously, ISSR and RAPD molecular markers were utilized to investigate intraspecific and intrageneric genetic variations in several Nepeta taxa [70][71][72][73][74].DNA barcoding regions, including nuclear (ITS) and plastid (e.g., matK, rbcL, trnL-F, trnS-trnG, and psbJ-petA) markers, were often utilized to reconstruct phylogenetic relations within the genus Nepeta [75][76][77].Petrova et al. [10] performed DNA barcoding by analyzing the N. nuda ITS sequences and regions of plastid DNA (rbcL/matK/trnH).Comparison of the obtained ITS sequences of N. nuda with the corresponding sequences of other Nepeta species revealed that it closely clustered with N. sheilae, N. deflersiana, N. isaurica, N. congesta, N. heliotropifolia, N. schiraziana, and N. cataria.However, three chloroplast markers showed its close relation to N. italica, N. tuberosa, N. cataria, N. grandiflora, and N. hemsleyana.Although both ITS and plastid markers have been proved to be suitable for determining taxonomic relationships among Nepeta species, they are considered too conservative to represent genetic variability within species and have thus not been considered for the present study.
The total number of alleles across the nine microsatellite markers used within the present study ranged from 2 for Cont039-gene0.9 to 12 for MN26.The basic parameters of genetic variation within and among the studied N. nuda populations are presented in Table 3.
Table 3. Localities of collection of Nepeta nuda material used in the study.The number of individuals per population used for the population genetics analysis and the yield of essential oils extracted from plants of the specific populations are also given (% PL-percentage of polymorphic loci; Namean number of alleles in a population; Ne-number of effective alleles; PA-number of private alleles; I-Shannon's informative index; Ho-observed heterozygosity; He-expected heterozygosity; uHe-unbiased expected heterozygosity; F-fixation index.).The mean number of alleles in a population ranged from 2.222 in the population Selište to 3.444 in the population Židilje, which is slightly higher than in Origanum compactum L. from Morocco [78] but considerably lower than in Rosmarinus officinalis L. from eastern Iberian Peninsula [79] and Salvia officinalis L. from the Balkans [80].The number of private alleles was low and ranged from 0 to 3, the latter being scored in the population Vlasina, which is in the range of the results obtained for several species from the same family [78,80].Observed heterozygosity ranged from 0.395 in Debeli Lug to 0.544 in Židilje, higher than that reported for O. compactum [57] but considerably lower than for R. officinalis [79] and S. officinalis [80].This parameter revealed a high heterozygosity level within individuals.Expected heterozygosity was generally slightly lower, with a span from 0.358 (uHe = 0.386) in Brodica to 0.524 (uHe = 0.589) in Gornje Selo, while Shannon's information index varied from 0.593 in Brodica to 0.931 in Gornje Selo.All the populations showed heterozygote excess, having negative F values, and only Debeli Lug had a slight heterozygote deficiency (F = 0.05).

Population
The AMOVA across all loci revealed that the greatest amount of genetic variance is captured within populations (80%), and among-population variance contributed 20% to the total variation (Figure 5A).This is quite an ordinary rule for outcrossing species belonging to the Lamiaceae family, as 78% and 84.73% within-population variation was recorded for O. compactum [78] and S. officinalis [80], respectively.Weir and Cockerham's [81] estimation of F ST was 0.203, which was significantly different from zero and similar to the value generated by AMOVA.
Pairwise F ST values between populations are presented in Figure 5B in colors ranging from green (no genetic differentiation) to orange (high genetic differentiation).Each value with F ST > 0.15 may imply significant population differentiation [64].It can be seen that the populations Selište and Donji Krivodol are rather strongly differentiated from the remaining populations, although they were not significantly different based on the chemical composition.
In Figure 5C, a two-dimensional PCoA scatterplot displays a grouping of individuals according to the population affiliation.It can be noted that the populations Vlasina, Gornje Selo, Brodica, and Debeli Lug failed to group along the two coordinates (Figure 5C).
The Mantel test revealed weak and a non-significant correlation (r = 0.3019, p ≤ 0.119) between geographic and genetic distances.STRUCTURE Harvester illustrated that K = 4, K = 6, and K = 7 are the most likely scenarios.In Figure 5D, the scenario with four genetic pools is presented.It can be seen that five out of six individuals from the population Selište draw mostly from a separate genetic cluster (blue), and that all individuals from Janjska reka almost completely consist of only one genetic cluster (green).The dominant genetic cluster in most of the remaining individuals is that marked red, while all the individuals from the populations Donji Krivodol and several individuals from Brodica, Rtanj, Židilje, Gornje Selo, and Vlasina draw mostly from the genetic pool marked yellow.
The analyzed set of EST-SSR and SSR markers proved not only to be very useful for the genetic diversity studies at the intra-species level but could possibly be transferred to other taxa of the genus Nepeta.The developed markers can further assist the germplasm conservation of N. nuda across its area.

Collection and Preparation of Plant Material
Aerial parts of flowering N. nuda plants were collected in June-July 2022 from 11 populations across Serbia (Figure 1).Plants were identified in the field by the authors.One plant from each population (except for the populations Rtanj, Židilje, and Gornje Selo) was stored in the Herbarium of the Institute of Botany and Botanical Garden, University of Belgrade, BEOU (acronym follows [82] and voucher numbers are presented in Table 4, along with the corresponding locations' names and their geographic coordinates).For DNA extraction and the preparation of methanol, extracts 2-3 leaves from at least 5 (mostly 10) plants per population were separately disposed into plastic zip-bags containing silica gel.Thus, each population was represented by 5 to 10 biological replicates (individuals/genotypes).Samples were kept in the dark at room temperature until use.For the isolation of essential oils, aerial parts (leaves, stems, and flowers) belonging to each population were pooled, air-dried at room temperature until constant weight was achieved, and subsequently mechanically chopped.Acetonitrile (Fisher Scientific, Leicestershire, UK) and formic acid (Merck, Darmstadt, Germany) were of MS grade.Ultra-pure deionized water was generated using the Water Purification System (New Human Power I Integrate, Human Corporation, Republic of Korea).Standards of cis,trans-nepetalactone, trans,cis-nepetalactone, 1,5,9-epi-deoxyloganic acid, and 5,9-dehydronepetalactone were isolated from natural sources as previously described by [19].Nepetanudoside A was quantified based on the calibration curve of loganin.Commercially available analytical standards of loganin, rosmarinic acid, caffeic acid, chlorogenic acid, and rutin (Sigma Aldrich, Hamburg, Germany) were used for quantification purposes.

Preparation of N. nuda Methanol Extracts
Silica gel-dried leaves (100 mg) of N. nuda were ground in liquid nitrogen and extracted with 1 mL of 96% methanol (w:v = 1:10).The extraction was performed by vortexing for 1 min and subsequent ultrasound-assisted extraction for 1 h at 4 • C. Supernatants, separated after centrifugation at 10,000× g for 10 min, were filtered through cellulose filters with 0.2 µm pore size (Agilent Technologies, Santa Clara, CA, USA).

Characterization of N. nuda Leaf Metabolites Using UHPLC-ESI-QToF-MS Analysis
One randomly selected N. nuda sample from each population was subjected to metabolic fingerprinting adopting a high-resolution mass spectrometry approach (HRMS).Analyses were performed using an Agilent 1290 Infinity UHPLC system in combination with a quadrupole time-of-flight mass spectrometer (6530C Q-QToF-MS, Agilent Technologies, Inc., Santa Clara, CA, USA).The mass detector was equipped with an electrospray ionization (ESI) source that was operated in both positive and negative ionization modes in the mass range from 100 to 1000 m/z.Chromatographic separation was performed at 40 • C on a Zorbax C18 column (2.1 × 50 mm, particle size 1.8 µm; Agilent Technologies, Inc., Santa Clara, CA, USA).The gradient elution program, UHPLC and MS parameters, as well as ion source settings, were described by Kostić et al. [83] and implemented in analyses of Nepeta species [66].
Agilent MassHunter software was used for acquisition, control, and MS data collection from the 6530C QToF-MS instrument.For the evaluation and presentation of MS data, R Studio software (enviPick and xcms R packages) was used as previously reported [84].Metabolites were identified based on their monoisotopic masses and MS 2 fragmentation and confirmed using previously reported metabolomic data about Nepeta species [19,27,30,32].Accurate masses of components and fragment ions were calculated using ChemDraw software (version 12.0, CambridgeSoft, Cambridge, MA, USA).The CAS SciFinder-n database was used to search for chemical compounds by formulae and structures (https://scifinder-n.cas.org/(accessed on 12 January 2024)).

UHPLC/DAD/(±)HESI-MS 2 Targeted Metabolic Profiling
Targeted UHPLC/DAD/(±)HESI-MS 2 metabolomics approach was adopted to quantify individual phenolics and iridoids in N. nuda methanol extracts.Analyses were performed on a Dionex UltiMate 3000 UHPLC system coupled with a DAD detector and configured with a triple quadrupole mass spectrometer (TSQ Quantum Access Max, Thermo Fisher Scientific, Basel, Switzerland).Samples were chromatographically separated on a Hypersil gold C18 analytical column (50 × 2.1 mm, 1.9 µm particle size; Thermo Fisher Scientific, Waltham, MA, USA) using the mobile phase gradient and a flow rate previously described by [19].The mobile phase consisted of (A) water + 0.1% formic acid and (B) acetonitrile + 0.1% formic acid.The injection volume was 10 µL.The selected reaction monitoring (SRM) mode of the instrument was used for the quantification of the targeted compounds by direct comparison with the commercial standards.Calibration curves of pure standards revealed good linearity, with r 2 values exceeding 0.99 (peak areas vs. concentration).Method validation parameters are presented in Supplementary Table S3.The total amount of each phenolic and iridoid compound was evaluated by the calculation of its peak area and is expressed as µg per 100 mg leaf dry weight (DW).Essential oils were prepared from aboveground flowering plants via hydrodistillation in a Clevenger-type apparatus connected to a 2-l borosilicate glass flat bottom flask (Isolab, Eschau, Germany), as previously described by [85].Hydrodistillation was performed for 2 h, using 50-100 g plant material and 1.5 l deionized water (Water Purification System, New Human Power I, Human Corporation, Seoul, Republic of Korea).The essential oil extraction procedure yielded from 2.19 to 5.97 µL essential oil per g dry weight (DW), depending on the N. nuda population (Table 2).

GC-MS Non-Targeted Metabolomics of N. nuda EOs
Profiling of volatile compounds contained in EOs of N. nuda plants was performed using an Agilent 8890 gas chromatography (GC) system coupled with a Mass Selective Detector (5977B GC/MSD, Agilent Technologies, Santa Clara, CA, USA) and connected to an automated sample extraction and enrichment platform (Centri ® , Markes International Ltd., Bridgend, UK).Chromatographic separations were performed for 21 min on an HP-5MS column (30 m × 0.25 mm, 0.25 µm film thickness; Agilent Technologies, Santa Clara, USA), using helium (99.999%,The Linde Group, Ireland) as a carrier gas at a flow rate of 1.6 mL min −1 .The settings of the MS procedure were previously described by Aničić et al. [66].In short, the transfer line was heated at 280 • C and the detector temperature was set to 270 • C. Mass spectra were acquired in the positive EI mode (+70 eV), with the temperature of the EI source set to 280 • C. Column temperature was linearly programmed from 40 to 300 • C at a rate of 20 • C min −1 and held isothermally at 240 • C for the subsequent 5 min.Centri ® was operating in a direct HS mode.Following dilution of EOs in 96% methanol (5 µL ml −1 ), 10 µL of an EO solution was transferred into a headspace (HS) vial and subsequently incubated at 40 • C for 10 min with agitation at 300 rpm.Sampling was performed in a split mode (20:1), with split flow of 24 mL min −1 , and the injector temperature set to 250 • C. Analyses were performed in the SCAN mode, tracking the compounds within the range of 45 to 500 amu.The constituents of the reaction mixtures were identified by comparison of their mass spectra and retention times with those of the respective standards, as well as by comparison with the NIST05 library.

Statistical Analysis of Metabolomics Data
For hierarchical cluster analysis (HCA), the input variables were scaled from 0 to 1. Principal component analysis (PCA) was constructed using the Past 4 software (version 4.14; [25,86].HCA was performed based on either Spearman's or Pearson's method of cluster agglomeration, using the Morpheus software (https://software.broadinstitute.org/morpheus (accessed 12 January 2024)).Quantitative metabolomics data were subjected to Tukey's post hoc test (p < 0.05) of one-way ANOVA.

EST-SSR Mining
For the RNA-Seq analysis of N. nuda leaf trichomes, leaf samples were collected in June 2020 from Balta Berilovac (Stara planina, East Serbia).Trichomes were harvested from fresh leaves using the dry ice abrasion technique, as detailed by [87].Subsequently, total RNA was isolated from detached trichome samples, and libraries were prepared by Macrogen Europe BV (Amsterdam, The Netherlands) using TruSeq Stranded Total RNA LT Sample Prep Kit (Plant) on the Illumina NovaSeq 6000 platform.RNA-Seq data were searched for microsatellite motifs (SSRs) using the Imperfect SSR Finder software [88].Motifs of three-to-five nucleotides in size and a minimum of five contiguous repeat units were screened.Primer pairs were designed using Primer3Plus with default settings [89].

DNA Extraction and Optimization of EST-SSR Markers' Amplification
DNA was isolated from silica gel-dried leaf samples using a modified CTAB method [90].DNA concentration and purity were assessed spectrophotometrically by measuring the absorbance at 230, 260, and 280 nm (NanoPhotometer N60, Implen, München, Germany).
We selected 20 putative EST-SSR markers to test for their variation in N. nuda.In total, thirty-three individuals from eleven N. nuda populations (three individuals per population) were used as a panel for the validation of the developed EST-SSR markers.Optimization of primers' annealing temperatures and amplification of EST-SSR markers were carried out using Mastercycler Nexus Gradient thermal cycler (Eppendorf AG, Hamburg, Germany).PCR reactions were performed in a final volume of 10 µL, each reaction containing 30 ng template DNA, 5.2 µL DreamTaq Green PCR Master Mix (2×; Thermo Fisher Scientific, Karlsruhe, Germany), 0.2 µM forward and reverse primer each, BSA 0.6 µg µL −1 (albumin, bovine serum; Sigma-Aldrich, Hamburg, Germany), and additional MgCl 2 to a final concentration of 4.5 Mm (Thermo Fisher Scientific, Karlsruhe, Germany).All EST-SSR loci were amplified using the following PCR reaction program: 94 Twenty PCR-amplified EST-SSR loci of the 33 selected N. nuda individuals were run on 1% agarose gels in 1 × Tris/borate/EDTA (TBE) buffer at 1.5 V cm −1 for 1.5 h, stained with ethidium bromide, and visualized in an UV transilluminator (ST4 3026-WL/26M, Vilber Lourmat, Collégien, France) along with a 50 bp DNA Ladder (Thermo Fisher Scientific, Karlsruhe, Germany).Amplification of 7 out of 20 EST-SSR markers yielded fragments of the expected length (or longer) across the tested plants.These seven EST-SSR loci were used for further population genetic studies of N. nuda (Table 5).  1 Loci multiplexed in the first run of fragment analysis; 2 loci multiplexed in the second run of fragment analysis.

Genomic Microsatellite Markers
To broaden the SSR genotyping, we tested 12 microsatellite markers reported by [91] for Mentha spp. that reportedly amplified DNA of Nepeta sp.After careful analysis of the amplification success and reproducibility, as well as of variation in alleles, two SSR markers were selected for further analysis: Cont028-gene0.2 and Cont039-gene0.9 (Table 5).

Amplification of the Selected Microsatellite Loci, Fragment Analysis, and Data Processing
The final PCR amplification of 9 selected molecular markers (7 EST-SSRs and 2 SSRs) was performed using 0.2 µM forward and reverse primers (Invitrogen, UK), the former being 5 ′ -labeled with one of the four fluorescent dyes (Table 5).The PCR reaction program was the same as for the amplification optimization.PCR reactions were performed without multiplexing.
Two runs per individual were prepared for fragment analysis, each containing 4 or 5 markers (Table 5).Fragment analysis was performed on an ABI3730xl DNA Analyzer (Thermo Fisher Scientific, Foster City, CA, USA) using GeneScan 500 LIZ size standard.Data scoring was performed using Peak Scanner (Thermo Fisher Scientific, Foster City, CA, USA).The following parameters were determined for each population: percentage of polymorphic loci, the number of different alleles, the number of effective alleles, the number of private alleles, Shannon's informative index, observed heterozygosity, expected heterozygosity, unbiased expected heterozygosity, and fixation index.To determine the sources of genetic variation within and between populations, AMOVA was performed using Wright's F statistics [81,92].Pairwise FST values between populations were also calculated [81].Visualization of populations' grouping was made through PCoA, using codominant genotypic distances among individuals.Population genetic analyses, AMOVA, and PCoA plotting were made by employing the GenAlEx 6.5 platform [93].
A mantel test between the genetic and geographic distances among the populations was performed using the IBD (Isolation By Distance) software [94] with 1000 randomizations.Population structuring was assessed using STRUCTURE 2.3.4 [95,96] with a burn-in length of 10,000 and a Markov Chain Monte Carlo (MCMC) of 100,000 randomizations, while the STRUCTURE HARVESTER online platform [97], which calculates the highest value of the second-order rate of change (DeltaK) according to the method of [98], was used to detect the number of K groups that best fits the dataset.

Conclusions
Based on previous studies dealing with the chemical characterization of different N. nuda accessions, it can be concluded that high intraspecific chemodiversity exists in this taxon.One must bear in mind that such implications might arise from different environmental conditions, growth stages of plants, sampling and extraction procedures, analytical techniques, and groups of metabolites analyzed.Literature data contain no comprehensive studies dealing with either inter-or intra-population variability of N. nuda aimed at determining sources of chemical diversity and providing explanations of this plasticity from the perspective of eco-evolutionary dynamics.The present study is the first attempt to explain variation in both metabolite profiles and genetic background within eleven populations of N. nuda from the Central Balkans.The low level of inter-population chemodiversity implies that differential environmental conditions might be less important for shaping the metabolomes of N. nuda populations, while their genetic backgrounds are of essential significance.This, at least, is true for the limited geographical area covered by the present study; this picture might vary if the whole species' area was studied.
As well as improving our comprehension of both the genetic and the chemical diversity of indigenous populations of this species, within-and among-population variability was assessed to support and to facilitate germplasm utilization.The recorded multidimensional genetic and chemical diversity of N. nuda represents a source of variance, which was harnessed in the quest to develop high-resolution chemical markers for a massspectrometry-guided selection of productive genotypes.Screening the distribution of well-defined bioactive compounds belonging to the groups of phenolics (e.g., rosmarinic acid) and iridoids (both aglycones and glycosylated forms), within and among N. nuda populations, provides a tool for the selection of elite representatives of practical importance.Furthermore, the developed chemical markers may be of importance when exploring intra-individual chemodiversity, which may further help in deciphering the tissue-and organ-specific, as well as the developmentally regulated, patterns of changes in metabolic fingerprints of Nepeta taxa, as influenced by endogenous (e.g., phytohormones) and exogenous (environmental) factors.The reconstruction of metabolic networks and the adequate definition of the role of different metabolic classes in the defense strategy of plants, as well as the explanation of their mutual effects, will inform out background knowledge of the factors that shape the metabolomes of plants and thus of their adaptive success in changing environments.Finally, the gathered knowledge on N. nuda chemical variation within the Central Balkans will help us to explain the molecular background of this diversity and the regulatory mechanisms involved, which is the course of our future work.

Figure 1 .
Figure 1.Map representing the populations of Nepeta nuda L. originating from the central Balkan Peninsula analyzed in the present study.For the coordinates, please refer toTable 1. Plant photos denote N. nuda individuals captured on site.

Figure 1 .
Figure 1.Map representing the populations of Nepeta nuda L. originating from the central Balkan Peninsula analyzed in the present study.For the coordinates, please refer to Table 4. Plant photos denote N. nuda individuals captured on site.

Figure 2 .
Figure 2. (A) Chemical structures of characteristic iridoid glycosides identified in Nepeta nuda methanol extracts (Glc-glucose or another hexose).(B) Heatmap of the scaled QToF MS data with the samples arranged according to the hierarchical cluster analysis (Spearman s method of cluster agglomeration).(C) Principal component analysis (PCA) biplot constructed based on the QToF MS data, with the first two PCs explaining 49.04% of the total variance among N. nuda accessions.Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots.

Figure 2 .
Figure 2. (A) Chemical structures of characteristic iridoid glycosides identified in Nepeta nuda methanol extracts (Glc-glucose or another hexose).(B) Heatmap of the scaled QToF MS data with the samples arranged according to the hierarchical cluster analysis (Spearman's method of cluster agglomeration).(C) Principal component analysis (PCA) biplot constructed based on the QToF MS data, with the first two PCs explaining 49.04% of the total variance among N. nuda accessions.Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots.

Figure 3 .
Figure 3. (A) Representative TIC GC/MS chromatograms of 11 samples of EOs originating from Central Balkans populations of Nepeta nuda.(B) Heatmap of the scaled GC/MS data with the samples (both populations and compounds) arranged according to the hierarchical cluster analysis and adopting the Spearman's method of cluster agglomeration.(C) Principal component analysis (PCA) biplot constructed based on the GC/MS data, with the first two PCs explaining 84.48% of the total variance among N. nuda accessions.Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots.

Figure 4 .
Figure 4. (A) UHPLC/(±)HESI-MS 2 quantitative data of targeted iridoids (green bars) and phenolics (orange bars) in methanol extracts of Nepeta nuda originating from 11 Central Balkans populations.Mean values ± SE are presented.(B) Principal component analysis (PCA) biplot constructed based on the quantitative data, with the first two PCs explaining 99.52% of the total variance among N. nuda accessions.For eleven N. nuda populations, 95% confidence ellipses are presented, colored according to the symbol legend.Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots.(C) Heatmap of the scaled UHPLC/(±)HESI-MS 2 quantitative data with populations and compounds arranged according to the hierarchical cluster analysis and adopting the Pearson s method of cluster agglomeration.

Figure 4 .
Figure 4. (A) UHPLC/(±)HESI-MS 2 quantitative data of targeted iridoids (green bars) and phenolics (orange bars) in methanol extracts of Nepeta nuda originating from 11 Central Balkans populations.Mean values ± SE are presented.Values labeled with different letters are significantly different (p < 0.05) according to post hoc Tukey's test of one way ANOVA.(B) Principal component analysis (PCA) biplot constructed based on the quantitative data, with the first two PCs explaining 99.52% of the total variance among N. nuda accessions.For eleven N. nuda populations, 95% confidence ellipses are presented, colored according to the symbol legend.Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots.(C) Heatmap of the scaled UHPLC/(±)HESI-MS 2 quantitative data with populations and compounds arranged according to the hierarchical cluster analysis and adopting the Pearson's method of cluster agglomeration.

Figure 5 .
Figure 5. Genetic diversity of Nepeta nuda populations from the Central Balkans, as estimated by EST-SSR and SSR markers.(A) Analysis of molecular variance.(B) Pairwise FST values of eleven populations of N. nuda.The green-to-orange color scale denotes low-to-high genetic differentiation.(C) PCoA biplot representing genetic clustering of individuals belonging to eleven populations of N. nuda.Cumulative variation explained by the two coordinates is 29.33% (19.34%COORD.1 + 9.99% COORD.2). (D) Individuals recruitment of the four genetic clusters inferred by STRUCTURE.

Figure 5 .
Figure 5. Genetic diversity of Nepeta nuda populations from the Central Balkans, as estimated by EST-SSR and SSR markers.(A) Analysis of molecular variance.(B) Pairwise F ST values of eleven populations of N. nuda.The green-to-orange color scale denotes low-to-high genetic differentiation.(C) PCoA biplot representing genetic clustering of individuals belonging to eleven populations of N. nuda.Cumulative variation explained by the two coordinates is 29.33% (19.34%COORD.1 + 9.99% COORD.2). (D) Individuals' recruitment of the four genetic clusters inferred by STRUCTURE.

Table 1 .
UHPLC-ESI-QTOF-MS data on metabolites identified in leaf extracts of eleven Nepeta nuda populations.
a Formic acid adduct; t R -retention time (min); ∆ mDa-mean mass accuracy; NA-not available.The [M ± H] ± column indicates in which ionization mode the corresponding compound was identified.Dominant fragments are represented boldfaced in the column MS 2 fragments.

Table 2 .
GC/MS characterization of aboveground extracts of 11 Nepeta nuda populations from the Central Balkans (relative content, %).

Table 4 .
Localities of collection of Nepeta nuda material used in the study.The number of individuals per population used for the population genetics analysis and the yield of essential oils extracted from plants of the specific populations are also given.
• C for 10 min; 36 cycles of 94 • C for 1 min; 52.5-61.0• C (depending on the calculated annealing temperature for primer pair, see Table 2) for 45 s; 72 • C for 45 s; and 72 • C for 10 min as a final extension step.

Table 5 .
Characteristics of nine microsatellite markers used for assessing the population structure of eleven populations of N. nuda.The expected sizes of the fragments (in bp) and their size ranges are given for a sample of 92 individuals.Ta = annealing temperature.