Variation in the chemical profiles of three foxglove species in the central Balkans

The aim of this study was to determine intra- and interspecies variation in the qualitative and quantitative composition of methanol-soluble metabolites in the leaves of three Digitalis species (D. lanata, D. ferruginea, and D. grandiflora) from the central Balkans. Despite the steady use of foxglove constituents for human health as valuable medicinal products, populations of the genus Digitalis (Plantaginaceae) have been poorly investigated to describe their genetic and phenetic variation. Following untargeted profiling using UHPLC-LTQ Orbitrap MS, by which we identified a total of 115 compounds, 16 compounds were quantified using the UHPLC(–)HESI–QqQ-MS/MS approach. In total, 55 steroid compounds, 15 phenylethanoid glycosides, 27 flavonoids, and 14 phenolic acid derivatives were identified across the samples with D. lanata and D. ferruginea showing a great similarity, while 15 compounds were characteristic only for D. grandiflora. The phytochemical composition of methanol extracts, considered here as complex phenotypes, are further examined along multiple levels of biological organization (intra- and interpopulation) and subsequently subjected to chemometric data analysis. The quantitative composition of the selected set of 16 chemomarkers belonging to the classes of cardenolides (3 compounds) and phenolics (13 compounds) pointed to considerable differences between the taxa studied. D. grandiflora and D. ferruginea were found to be richer in phenolics as compared to cardenolides, which otherwise predominate in D. lanata over other compounds. PCA revealed lanatoside C, deslanoside, hispidulin, and p-coumaric acid to be the main compounds contributing to the differences between D. lanata on one side and D. grandiflora and D. ferruginea on the other, while p-coumaric acid, hispidulin, and digoxin contribute to the diversification between D. grandiflora and D. ferruginea. However, quantitative variation in the metabolite content within species was faint with mild population diversification visible in D. grandiflora and particularly in D. ferruginea. This pointed to the highly conserved content and ratio of targeted compounds within the analyzed species, which was not severely influenced by the geographic origin or environmental conditions. The presented metabolomics approach might have, along with morphometrics and molecular genetics studies, a high information value for further elucidation of the relationships among taxa within the genus Digitalis.


Introduction
Both genotyping and phenotyping of a robust sample set of a species under study are essential to properly understand the variation in its populations. The variation itself is the main prerequisite for a species to "outwit" environmental pressures and to prolong its survival on the planet. In medicinal plants, assessing the genetic and phenetic variation of wild plant populations enables appropriate planning of selective breeding to ensure favorable yield. Besides crop and timber improvement programs, medicinal plants, being a great part of the modern world health care system (Fabricant and Farnsworth, 2001), are in the focus of plant breeders. In fact, around 25% of drugs commonly used worldwide is derived from plants (Mishra et al., 2013).
During the Middle Ages, foxglove (Digitalis spp., family Plantaginaceae) has been considered to have multifarious healing properties (Aronson and Aronson, 1986) and the first reports of its role as an agent in the congestive heart failure treatment dated from the year 1250 (Williams, 1861). Today, D. purpurea L. and, more notably, D. lanata Ehrh. are both widely collected from nature and field-cultivated to extract cardenolides, being among the best producers of these compounds in the plant world (Clemente et al., 2011). One of them, digoxin, is mainly produced in Europe from dried leaves of D. lanata, reaching up to 1.5% dry weight (Kennedy, 1978;Bown, 1995;Clemente et al., 2011). There is a vast diversity of cardenolides within the genus Digitalis with more than one hundred different compounds isolated, the most commercially attractive, besides digoxin, being lanatoside C, digitoxin, and acetyldigitoxin (Clemente et al., 2011;Kreis, 2017). Other compounds, such as digitanols, steroidal saponins, anthranoids, phenols, sterols, and polysaccharides, have been identified in various Digitalis species (see Clemente et al., 2011 andKreis, 2017 for reviews). Today, with the development of synthetic substitute drugs for cardiac diseases, foxgloves are given "a second chance" (Kreis, 2017) as potent agents against various viruses as well as for the treatment of cystic fibrosis and several cancer types (Bertol et al., 2011).
Out of the 27 recognized species of the genus Digitalis (Plants of the World Online | Kew Science, 2023), 5 are found in the Balkan Peninsula: D. lanata, D. laevigata Waldst. & Kit., D. ferruginea L., D. viridiflora Lindl., and D. grandiflora Mill. (syn. D. ambigua Murray), the first three belonging to the section Globiflorae and the last two to expanded "Maranthae" (Bräuchler et al., 2004). Two more species, namely D. ikarica (P.H.Davis) Strid and D. fuscescens Waldst. & Kit., can be found in the Balkans, the former inhabiting several North Aegean and Dodecanese islands, while the latter being reported to represent a hybrid between D. grandiflora and D. purpurea (Braemer et al., 1927). Among these, D. lanata is somewhat better studied in terms of genetic and phenetic diversity and by its potential to accumulate various specialized metabolites (Braga et al., 1997;Yücesan et al., 2018). Few data are available for the other species listed (Kennedy, 1978;Katanićet al., 2017;Kaska et al., 2020). Surprisingly, variation in metabolite fingerprints both between and within wild populations is pretty much neglected, having in mind that it can render starting material for selective breeding or at least for the selection of specific genotypes that would accumulate metabolite amounts appropriate for the pharmaceutical industry. Nevertheless, according to EURISCO (Kotni et al., 2023), many landraces are stored in seed collections. Therefore, phytochemical characterization of an appropriate sample set consisting of as many populations as possible and represented by a sufficient number of individuals is an evident prerequisite to perceive the extent of variation of specialized metabolites in Digitalis species. Since the Balkan Peninsula represents one of the two centers of the genus diversity and considering that D. lanata, D. ferruginea, and D. grandiflora are most commonly found throughout it, the three species were selected for the study. In a greater perception, D. lanata grows across the southeast Europe and Turkey (and has been introduced into North America and several central European and Asian countries), D. ferruginea, besides the Balkan Peninsula, inhabits the Apennine Peninsula, Asia Minor, and the Caucasian region, while D. grandiflora is spread throughout Europe (except for several westernmost and northernmost countries) and northwest Asia, but has been introduced into North America. Central Balkan populations of these three species contain few individuals (several dozen at maximum, pers. obs.), which implicates that they might be at risk of extinction due to the bottleneck effect. Therefore, the protection of degraded populations should be considered and the level of their genetic variation must be assessed (Nebauer et al., 2000). However, only D. grandiflora has been studied in this sense (Boronnikova et al., 2007). The detection of population differentiation by studying either genetic or phenetic variation (or preferably both) provides valuable information for planning conservation measures and conducting monitoring as well as for sampling germplasm for ex situ conservation (Nazir et al., 2008;Clemente et al., 2011).
The present study is aimed at determining, disentangling, and interpreting the chemical diversity of three Digitalis species growing in the Balkan Peninsula, to setup the footprint for meaningful biodiversity conservation strategies and to propose the means for sustainable utilization of bioresources. Furthermore, this study will establish the background for future elucidation of molecular aspects of chemical diversity by adopting state-of-the art omics technologies and will enable the development of alternative strategies for the production of bioactive compounds through biotechnology approaches. The specific objective of this study was to select high-resolution chemical markers to estimate both interand intrapopulation variability of Digitalis species by applying a targeted metabolomic approach.

Chemicals
All reagents and solvents used were of analytical grade. Acetonitrile, formic acid (both MS grade), and methanol (HPLC grade) were purchased from Merck (Darmstadt, Germany). Ultrapure water (Water Purification System, New Human Power I Integrate, Human Corporation, Republic of Korea) was used to prepare standard solutions and blanks. Analytical standards of protocatechuic acid, syringic acid, p-hydroxybenzoic acid, 5-Ocaffeoylquinic acid, caffeic acid, aesculetin, isoorientin, p-coumaric acid, quercetin 3-O-glucoside, naringin, luteolin, hispidulin, and isorhamnetin were purchased from Sigma Aldrich (Steinheim, Germany). Standards of deslanoside, lanatoside C, and digoxin were provided by Professor Yang Ye and Dr. Chunping Tang (Shanghai Institute of Materia Medica-SIMM, Chinese Academy of Sciences -CAS, China).

Plant material
Leaves from the flowering stems were collected from individual plants during June, July, and August of 2020 and 2021 from wild populations in the central part of the Balkan Peninsula. Localities and the list of 28 populations of three Digitalis species (Digitalis grandiflora, D. lanata, and D. ferruginea) are listed in Table 1 and their spatial distribution is presented in Figure 1. Plants were identified in the field by the authors and classified by referring to The World Flora Online (WFO) database (The World Flora Online, 2023)

Methanol extracts preparation
After harvesting, leaves were immediately transferred into plastic zip-bags containing silica gel. Prior to methanol extraction, leaves were ground to a fine powder using liquid nitrogen. Approximately 50 mg of dry plant material was extracted with 1 ml of 80% methanol overnight at room temperature. The next day, samples were sonicated (Sonorex Bandelin Electronic, Berlin, Germany) for 1 h and subsequently centrifuged at 10,000g for 10 min. The supernatants were filtered through 0.2-mm cellulose filters (Agilent Technologies, Santa Clara USA) and stored at 4°C until use.

Identification and quantification of metabolites 2.4.1 UHPLC-LTQ Orbitrap MS untargeted metabolomics analysis
Identification of metabolites in methanol extracts of the three foxglove species was done by an untargeted approach using an Accela UHPLC system connected to a linear ion trap-Orbitrap hybrid mass spectrometer (LTQ OrbiTrap XL, Thermo Fisher Scientific, Bremen, Germany) with heated electrospray ionization (HESI). Methanol extracts of D. lanata, D. ferruginea, and D. grandiflora were prepared from leaves collected in 2020 in Zavojsko jezero (population accession code No. DLZAV), Zaovine (population accession code No. DFZAO), and Debelo brdo (population accession code No. DADB), respectively (Table 1).
The mass spectrometer was operated in either negative or positive ionization mode, depending on the compound class. The HESI-source parameters were described by Koprivica et al. (2018). MS spectra were acquired by full range acquisition covering 100-1500 m/z. The data-dependent MS/MS events were always performed on the most intense ions detected in the full scan MS. The ions of interest were isolated in the ion trap with an isolation width of 5 ppm and activated with 35% collision energy levels.
Metabolites were identified according to the corresponding spectral characteristics: mass spectra, accurate mass, characteristic fragmentation patterns, and corresponding retention time. Tentative identification of various metabolites was achieved by studying their MS n spectra and comparing them with the available literature on spectroscopic and mass data for compounds detected previously in the genus Digitalis and other related species (Supplementary Table 1). Xcalibur software (version 2.1) was used for the instrument control, data acquisition, and data analysis.

UHPLC/(-)HESI-QqQ-MS/MS targeted metabolomics analysis
Quantification of the targeted compounds was performed using a Dionex Ultimate 3000 UHPLC system connected to a triplequadrupole (QqQ) mass spectrometer (TSQ Quantum Access Max, Thermo Fisher Scientific, Bremen, Germany). Metabolite quantification was performed across 259 individual plants belonging to the three study species (Table 1).
A Syncronis C18 analytical column (100 × 2.1 mm) with 1.7 µm particle size (Thermo Fisher Scientific, Bremen, Germany) was used for the chromatographic separation. The flow rate and the composition of the mobile phases as well as the gradient elution program are described in the previous section (2.4.1.). The mass detector was equipped with a HESI source operated in the negative ionization mode. The parameters of the HESI source and the other mass detector settings were previously described by Banjanac et al. (2017). The selected reaction monitoring (SRM) mode of the instrument was used for the quantification of the targeted compounds in the samples. The compounds were quantified by direct comparison with commercial standards. Calibration curves revealed good linearity, with r 2 values exceeding 0.99 (peak areas vs. concentration). The total amount of each compound was evaluated by calculation of the peak area and is expressed as mg/kg. The concentration of digoxin is expressed via the calibration curve of deslanoside, because the available amount of digoxin standard was insufficient to obtain the calibration curve with appropriate calibration levels.

Statistical analysis
The quantitative data were analyzed by two unsupervised methods, principal component analyses (PCA) and hierarchical cluster analyses (HCA), as well as one supervised learning method, linear discriminant analysis (LDA). For HCA the input variables were scaled to the [0, 1] range. HCA was performed based on Euclidean distances with cluster agglomeration using Ward's (Ward, 1963) minimum variance method. LDA was performed with a presumption of the equal prior probability of classes. The correlation matrix for the quantitative data was constructed using Pearson's correlation coefficients. All statistical analyses were performed in the Past 4 software (version 4.12; Hammer et al., 2001).

Results and discussion
Although some data on the phytochemical content of D. grandiflora, D. ferruginea, and especially D. lanata are available in the literature, these species are only scarcely investigated for the phytochemical interpopulation variability (Yücesan et al., 2018). Several previous studies have reported carbohydrates, iridoids, and caffeoyl phenylethanoid glycosides as suitable chemomarkers for the Plantaginaceae family, including the genus Digitalis (Taskova et al., 2005). Cardenolides and phenolics were not considered in this regard. On the other hand, tracing the amounts of commercially attractive compounds can provide valuable information to study the diversity of Digitalis species in the phylogenetic context. This may also highlight their characteristic populations' metabolite profiles allowing to pick out the genotypes that produce pharmaceutically interesting compounds. More than 100 different cardenolides have been isolated from the genus Digitalis, and previous studies have been focused mainly on several well-known and exploited Digitalis species (mostly D. purpurea and D. lanata). However, to search for novel compounds, there is a necessity for in-depth screening of various Digitalis species for their metabolite composition, which may render bioactive agents for the treatment of various diseases. On the other hand, cardenolides occur in multicompound mixtures whose composition represents complex phenotypes that may vary along multiple dimensions and levels of biological organization. Digitalis species are characterized by the presence of a vast array of phenolic compounds, among which flavonoids, mainly belonging to the flavone and 3-methoxyflavone groups, predominate. These two large groups of bioactive compounds and their diversity both within and among Digitalis species deserve a considerable attention. Simultaneous tracing of their composition in samples can ensure characterization of phenotypes in an ecologically and evolutionary meaningful way and enable the elucidation of the relations among and within group of compounds pointing out to the links between their biosynthetic pathways.

Species-specific metabolite profiles
According to Kreis (2017), HPLC-MS procedures for assessing the metabolic footprints of cardenolides found in Digitalis spp. are highly recommended. Identification of metabolites present in the methanol extracts of three foxglove species (D. lanata, D. ferruginea, and D. grandiflora) was achieved by high-resolution mass spectrometry (HRMS) in combination with MS n fragmentation. The representative base peak chromatograms of the three species are shown in Supplementary Figure 1. Using this technique, the molecular formula of an unknown compound can be determined through the exact mass, and its structure can be proposed and solved by studying its fragmentation pathway. Using the UHPLC-LTQ OrbiTrap MS technique in both positive and negative ionization modes, totally 115 compounds were identified based on their monoisotopic masses, MS n fragmentation, and previously reported MS data ( Map presenting the populations of Digitalis grandifloraa, D. lanata, and D. ferruginea originating from the central Balkan Peninsula analyzed within the present study. For the population labels, please refer to Table 1. flavonoid aglycones (9 compounds), phenolic acid derivatives (14 compounds), and 6 compounds belonging to other classes. In addition to the expected cardiac glycosides, steroidal saponins as well as pregnane and furostanol glycosides were also present in the analyzed extracts. A detailed LC/MS qualitative analysis revealed very similar profiles of steroidal glycosides in the extracts of D. lanata and D. ferruginea leaves, in particular of digoxin, deslanoside, and lanatosides A, B, and C, which is consistent to earlier reports (Kennedy, 1978;Braga et al., 1997;Verma et al., 2014;Yücesan et al., 2018). On the other hand, these compounds were not detected in leaves of D. grandiflora although digoxin has been reported in small amounts in an earlier study (Kennedy, 1978). Table 2 specifies which of the identified compounds were previously detected in any Digitalis species, while for a certain number of compounds it can be concluded that they were identified for the first time in Digitalis species.
The largest number of identified compounds (55 in total) belongs to the group of steroids (cardenolides, pregnane glycosides, furostane-type and spirostane-type steroidal saponins), 43 of which are steroidal glycosides and 12 are steroidal aglycones. Different sugar derivatives can be found as sugar components of these cardenolides and other steroid glycosides. The most common is digitoxose, which gives 130 Da as a neutral loss in the mass spectrum, followed by different hexoses (162 Da), acetyldigitoxose (172 Da), digitalose (160 Da), and less often deoxyhexose (146 Da) (Ravi et al., 2020). As can be seen from Table 2, some of these steroid derivatives are named after a certain Digitalis species (e.g., digitalin, purpureagitoside, lanatoside, purpurea glycoside), while some compounds are named trivially as furostanol or spirostanol glycosides. All derivatives of steroidal aglycones are, due to their polarity, identified only in the positive ionization mode. For all listed compounds (except compound 23) from the group of steroids, Supplementary Table 1   base peak at 311 m/z. Interestingly, a compound with this structure has not yet been identified but has been synthesized as a digitoxin-related compound and used to test its anticytomegalovirus activity (Cai et al., 2014). A total of 15 compounds was identified from the group of phenylethanoid glycosides (PG), all of which except decaffeoyl acteoside isomers (compounds 56 and 57) were previously identified in various Digitalis species (Table 2). All the three examined Digitalis species proved to be rich in PG, as previously reported (Brieger et al., 1995), except for compound 69 (forsythiaside), which was identified only in the extract of D. grandiflora. Compounds from the group of phenylethanoid glycosides specific for Digitalis have been shown to be potential tumor inhibitors -cytotoxic activity of the PG isolated from D. davisiana showed a strong impact on HEp-2 cells (Kutluay et al., 2019), while PG extracts from D. purpurea showed PKCainhibitory bioactivity (Zhou et al., 1998).
All 18 identified flavonoid glycosides belong to the flavone subgroup and are represented by hydroxyl, methyl, and methoxy derivatives of apigenin and luteolin (Table 2). Compound 71 (apigenin 6,8-di-C-hexoside) is a C-glycoside, while all others (compounds 72-88) are identified as 7-O-glycosides. Flavonoid C-glycosides were not found to be common for Digitalis species, but compound 71 was previously identified in the related genus Plantago, belonging to the same family Plantaginaceae (Kawashty et al., 1994). Among the derivatives of flavonoid Oglycosides, hexosides (neutral loss of 162 Da), hexuronides (neutral loss of 176 Da), and two acylhexosides (compound 83 and 88) were found. The presence of three glycosides (compounds 73, 76, and 77) with two hexuronic acid units was confirmed by the specific MS 2 fragmentation, in which the fragment 351 m/z (dihexuronyl-H) was formed as a base peak. Usually, the MS 2 base peak actually represents the mass of the deprotonated aglycone, which was not the case here, but the disaccharide residue had the highest abundance. In the next fragmentation stage, the MS 3 base peak at 193 m/z was formed by the loss of 158 Da (dihexuronyl-H 2 O). A similar fragmentation pathway was noticed for compound 73 (luteolin 7-O-hexosylhexuronide), with the difference that its MS 2 base peak corresponds to the mass of the deprotonated aglyconeluteolin (285 m/z). It was previously reported that leaves of D. purpurea contained luteolin 7-O -hexosyl-hexuronide (Harborne, 1963). Compound 86 (jaceosidin 7-O-hexuronide) has not been detected in Digitalis species so far, but its aglycone jaceosidin (Imre et al., 1984) and 7-O-glucoside of jaceosidin (Hiermann et al., 1977) have been reported. Using the exact mass of compound 88 with molecular ion at 563.14063 m/z in the positive ionization mode, its molecular formula C 26 H 26 O 14 was generated. By studying its fragmentation pattern, it was concluded that dihydroxy-dimethoxyflavone, malonic acid and hexose had to be present in its structure. In the first stage of fragmentation, the malonylhexose unit was lost and the MS 2 base peak at 315 m/z was formed, which corresponds to the mass of protonated dihydroxy-dimethoxyflavone. Further, in MS 3 fragmentation, the loss of the CH 3 group (15 Da) resulted in only one fragment at 300 m/z, while in MS 4 fragmentation, a base peak at 168 m/z was observed. The fragment ion at 168 m/z was formed by the specific RDA fragmentation of flavonoids and it can be designated as a 1,3 A + -CH 3 fragment (Tsimogiannis et al., 2007), which further leads to the conclusion that this flavonoid has two hydroxyl groups and one methoxy group on the A ring. In view of all the above-mentioned facts, 6-methoxy-4'-methylapigenin or pectolinaringenin was proposed as the aglycone of this compound, which was otherwise common for Digitalis species (Imre et al., 1984). A proposed structure and a fragmentation pathway of this compound, trivially assigned as pectolinaringenin 7-O-malonylhexoside, is depicted in Supplementary Figure 3.
Among the nine identified flavonoid aglycones (Table 2), eight were flavone derivatives and only one (compound 96) was identified as flavonolsantin. Santin belongs to the subgroup of flavonoids termed O-methylated flavonols and it was earlier isolated from D. orientalis (Imre et al., 1984). The presence of compounds 92-94 was confirmed by comparing their mass spectra with the analytical standards, while the other compounds from this group were proven by determination of the exact masses and fragmentation pathways (ApSimon et al., 1963;Hiermann et al., 1977;Imre et al., 1984).
Derivatives of phenolic acids found in this study (Table 2) were previously reported to be specific for Digitalis species and were actually derivatives of hydroxycinnamic acid (Kırmızıbekmez et al., 2009;Katanićet al., 2017;Kaska et al., 2020). In this study, a certain number of hydroxybenzoic acid derivatives (compounds 98, 100, and 102) were also found. Three compounds (102, 104, and 105) were identified by direct comparison with the corresponding standards. A significant number of hexosyl derivatives of phenolic acids were detected, giving a specific fragmentation with the neutral loss of 162 Da.
Six compounds that structurally do not belong to any of the above-mentioned groups are included in the "other compounds" category ( Table 2). The presence of two derivatives of coumarins, aesculetin (compound 110), and its 6-O-glucosideaesculin (compound 112) was confirmed by comparison with the analytical standards. Compound 111, hebitol II (O-acyl carbohydrate), was found only in the D. lanata extract in this study, and its fragmentation was consistent with previously published MS data (Frisčǐćet al., 2016). Compound 113, lignan syringaresinol, has not been previously detected in Digitalis species, but its fragmentation is in agreement with the literature data (Jiang et al., 2019). The last two compounds (114 and 115), digitoemodin and w-hydroxyziganein-1-methyl ether, are anthraquinones previously isolated from D. cariensis (Imre et al., 1994).

Quantitative patterns of metabolites found in three Digitalis species
Targeted metabolic profiling was conducted with the aim of quantifying 16 bioactive compounds in methanol extracts of the three analyzed Digitalis species. We were able to quantify 3 compounds from the class of cardenolides: digoxin, lanatoside C, and deslanoside (syn. desacetyllanatoside C). In addition, we quantified 7 phenolic compounds, derivatives of either hydroxycinnamic (5-O-caffeoylquinic acid, caffeic acid, p-coumaric acid, and aesculetin) or hydroxybenzoic acid (p-hydroxybenzoic acid, syringic acid, and protocatechuic acid). Three Digitalis species have also been investigated with respect to their flavonoid content by targeting totally 6 compounds belonging to flavones (luteolin, hispidulin, isoorientin), flavonols (quercetin, isorhamnetin), and flavanones (naringin). The results indicate that the most abundant phenolic acids of D. grandiflora belong to the group of hydroxybenzoic acids (protocatechuic, syringic, and phydroxybenzoic acid), while p-coumaric acid predominates among hydrohycinnamic acids (Gasǐćet al., 2023). Dominant phenolic acids in D. lanata and D. ferruginea are p-coumaric acid and phydroxybenzoic acid (Gasǐćet al., 2023), which slightly coincides with previous findings, reporting chlorogenic and caffeic acids (Katanićet al., 2017) or ferulic acid (Kaska et al., 2020) as the most abundant in D. ferruginea. Hispidulin and other flavonoids belonging to the flavone and 3-methoxyflavone groups are especially abundant in Digitalis species (Kaska et al., 2020), which is also confirmed in the present study. Naringin, luteolin, and hispidulin are the most abundant flavonoids in D. grandiflora (Gasǐćet al., 2023). In D. lanata hispidulin predominates, while D. ferruginea leaves are characterized by the prevalence of hispidulin and naringin.
The content of the main cardenolides in D. lanata varies with the growth stages (Freier, 1977;Braga et al., 1997) and is further influenced by the leaves' maturity (Vogel and Luckner, 1981), growth conditions, vegetation stage, collection time, drying method, and storage conditions and duration (Balakina et al., 2005). Within the present study, the profiling of specialized metabolites has been performed using silica geldried leaf samples collected from blooming plants in their second year of growth, as recommended by Clemente et al. (2011) at the stage when the content of pharmaceutically interesting compounds is at the optimal level. Leaves at the same developmental stage were collected to further reduce the effect of developmental stage on the content of compounds of interest. All the three targeted cardenolides were present in D. lanata, with lanatoside C predominating (Gasǐćet al., 2023). According to the literature, D. lanata is known to contain primary cardenolide glycosides including lanatoside A, B, and C series (Braga et al., 1997;Peŕez-Alonso et al., 2012;Yücesan et al., 2018). On the other hand, D. grandiflora has been reported to contain lanatoside A and B, but lanatoside C was absent (Wichtl et al., 1987). This is in accordance with the present study, since we were able to quantify digoxin in D. grandiflora leaves, while lanatoside C and deslanoside were not detected or were present in trace amounts (Gasǐćet al., 2023). Major cardenolides in leaves of D. ferruginea were digoxin and deslanoside, which coincides to a previous study (Verma et al., 2014) that reported digitoxigenin and digoxigenin as the dominant cardenolides in this species.
The content and the ratio of cardenolides in Digitalis species are regulated by the activity of enzymes responsible for the conversion of primary glycosides such as lanatoside A, lanatoside B, and lanatoside C into corresponding secondary glycosides digitoxin, gitoxin, and digoxin, respectively (Pellati et al., 2009). The hydrolysis of lanatoside C to form digoxin proceeds in two stages ( Figure 2A) and involves the enzyme digilanidase and one deacetylation step (Cobb, 1976). Deslanoside (desacetyllanatoside C) is a product of the lanatoside C metabolism and a precursor of digoxin ( Figure 2A). The hydrolysis of lanatoside C usually occurs after the damage of leaf tissue or plant harvesting. It has been suggested that the absence or low enzymatic activity in the plant material favors the preservation of primary glycosides in D. lanata (Balakina et al., 2005). On the contrary, if the enzymes are sufficiently active, a spontaneous enzymatic degradation of primary to secondary glycosides occurs, resulting in an extremely complex profile in which both primary and secondary glycosides are present in different ratios, depending on both genetic background and post-harvest processing methodologies (Balakina et al., 2005). Within the present study, harvested leaves were immediately stored on silica gel, which enabled their fast dehydration and prevented degradation of primary into secondary glycosides. Thus, we were able to trace the species-and genotype-specific content and ratio of the targeted cardenolides in D. lanata, D. ferruginea, and D. grandiflora leaves. Accordingly, we may conclude that D. lanata is characterized by the predominance of lanatoside C, which is moderately converted to digoxin, most probably via deslanoside (Figure 2A), as this intermediate is present in significant amounts in leaves of this species. The alternative digoxin precursor, acetyldigoxin, was not recorded in any analyzed taxon (Table 2). On the other hand, D. grandiflora leaves most likely display prominent activity of the enzymes responsible for the conversion of primary into secondary glycosides, as digoxin is found to be the major cardenolide in this species (Gasǐćet al., 2023). Lanatoside C and its metabolite deslanoside were absent in the majority of the analyzed D. grandiflora accessions or were present in trace amounts in several samples. In D. ferruginea, lanatoside C is, most likely, efficiently converted into digoxin, the dominant cardenolide compound. The conversion probably takes place via deslanoside intermediate, as this compound is the second most abundant compound in this species.

FIGURE 2
A simplified scheme of the lanatoside C hydrolysis that results in the formation of digoxin in two stages and involves the enzyme digilanidase and one deacetylation step (A). The conversion goes via either deslanoside (desacetyllanatoside C) or acetyldigoxin. D. lanata (a) is a rich source of lanatoside C, D. ferruginea (b) contains predominantly deslanoside and digoxin, while D. grandiflora (c) is abundant in digoxin. Correlation analysis of the compounds' quantitative data was performed (B), and the correlation matrix was constructed using Pearson's correlation algorithm. Color scale indicates a positive (blue) or a negative (red) correlation.  Table 1.

A B C
We further performed the pairwise correlation analysis by calculating the Pearson's correlation coefficient ( Figure 2B) to provide further evidence for the stated presumptions. Significant positive correlations were observed between the contents of lanatoside C and deslanoside. Digitoxin was negatively correlated with the other two cardenolides analyzed here, which is not surprising taking into account that this secondary glycoside emerges by the deacetylation of lanatoside C to deslanoside, which then undergoes hydrolysis. Cardenolides were further positively correlated with hydroxybenzoic acids and flavones hispidulin and isorhamnetin. Positive correlations were mainly observed within the group of phenolics, with the exception of syringic acid, which was negatively correlated with the majority of analyzed compounds.
To better understand the relationship between the analyzed taxa (D. grandiflora, D. lanata, and D. ferruginea), we performed PCA analysis on phytochemical dataset acquired for 16 targeted compounds ( Figure 3A). PC1 and PC2 cumulatively explained 83.56% of the total variance. D. lanata is clearly separated from the other two Digitalis species along the PC1, explaining 61.63% of the total variation. On the other hand, D. grandiflora and D. ferruginea segregate along the PC2, which explains 21.93% of the variability. The major contributors to PC1 are lanatoside C, deslanoside, hispidulin, and p-coumaric acid ( Figure 3B). Along the PC2, samples are distinguished mainly by p-coumaric acid, hispidulin, and digoxin ( Figure 3C). This unsupervised multivariate method offered a glimpse of possible usefulness of the selected combination of metabolites for phytochemical determination of the three Digitalis species. Although we presumed that D. lanata and D. ferruginea, phylogenetically closely related taxa (Bräuchler et al., 2004) would show greater similarity in phytochemical profiles, this was not the case.
Based on the HCA analysis (Supplementary Figure S4), a similar conclusion to that derived from PCA can be drawn regarding the phytochemical relationships among the analyzed Digitalis taxa. However, HCA more clearly depictured the populations' linkages and indicated clustering of populations based on their geographical distribution.

Intraspecies quantitative metabolite patterns
Having in mind that distinct metabolic pathways of the 16 compounds may be reflected on different accumulation patterns across the studied populations within the three Digitalis species, one would expect huge quantitative variation in their profiles, since interpopulation variation in metabolite accumulation has been previously reported for other Digitalis species (Nebauer et al., 1999;Usai et al., 2007). However, our results imply poor population differentiation in this regard for all the three studied species. As shown in Figures 4-6, a divergence of two populations from the central sample cloud in the PCA can be observed for D. ferruginea ( Figure 6A) and only mild diversification of several populations of D. grandiflora can be noticed in the LDA ( Figure 4B). In particular, all individuals representing the population DABARE and 5 individuals each from DAZ and DADZ cluster separately from the main sample cloud in the HCA matrix formed by the remaining samples belonging to 12 other populations of D. grandiflora ( Figure 4C). The three populations mentioned are from Dinaric karst areas with limestone as the basal rock, located at the high altitudes (Table 1), which could be the reason for their segregation from the common cluster, since soil characteristics are reportedly related to cardenolide production (Roca-Peŕez et al., 2004). The most abundant compound in the majority of D. grandiflora populations, p-hydroxybenzoic acid, is in the same time the compound that varies at the most, followed by two flavonoid aglycones, hispidulin and luteolin ( Figure 4A). Although present in greater amounts than in D. lanata, digoxin is evenly quantified across the populations of D. grandiflora (Gasǐćet al., 2023). As the species with the broadest areal among the three studied species, D. grandiflora shows rather surprisingly low variation in the metabolite profiles, despite being represented by the largest sample set. D. grandiflora is a biennial plant, occurring as a basal rosette of leaves during the first year, while the second-year rosette leaves wither rapidly during flowering stem elongation. To uniform the sampling procedure, leaves of the same developmental stage originating from flowering stems were collected for the analyses. Delayed flowering of D. grandiflora in higher altitudes enabled us to perform sampling throughout summer seasons and, presented results indicate that sampling period had no significant influence on the content of metabolites. An earlier study (Boronnikova et al., 2007) that estimated the partitioning of molecular variation showed that only a small fraction of the species' genetic variation was present between populations, while most of it resided within populations (about 90%, depending on the type of molecular marker used). Therefore, the weak population differentiation seen in the metabolite profiles would most certainly have its ground in their low genetic differentiation. Further studies on the genetic background of these populations currently planned by the same team of authors will undoubtedly clarify this hypothesis.
The sample matrix made of 6 populations of D. lanata presented in Figure 5A is much more homogeneous, while clustering fails to offer the resolution to distinguish potential groups. Hispidulin and p-coumaric acid are the major factors of the summary metabolite variation among the studied populations of this species ( Figure 5A). They are followed by the two most abundant cardenolide compounds in D. lanata, lanatoside C and deslanoside. The four compounds are also leading in the quantified amounts throughout the studied populations (Gasǐćet al., 2023). Nearly 50 years ago, D. lanata genotypes were actually quite often assessed but only with respect to cardenolide content (Weiler and Zenk, 1976;Weiler and Westekemper, 1979;Wichtl et al., 1987) and no Principal component analysis (PCA) biplot with the first two PCs explaining 84.38% of the total variance (A) among D. grandiflora accessions. For fifteen D. grandiflora populations (each labeled by a different symbol) 95% confidence ellipses are presented. Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots. Linear discriminant analysis (LDA) (B) represents the same fifteen populations of D. grandiflora with 95% confidence ellipses. Heatmap of the scaled quantitative data (C) with the samples arranged according to the hierarchical cluster analysis (Ward's method of cluster agglomeration). Compounds' labels are the same as in Figure 3. For the interpretation of population abbreviations in the figure legend, please refer to Table 1. comprehensive study based on a broader metabolite profiling has been performed in this species. However, these articles report on a remarkable variation in the accumulation of cardenolides in plants originating in different geographic regions. In our study, we were able to collect samples from 6 populations having the maximum distance of 200 km and growing in similar habitats, and only mild differences in metabolic profiles were anticipated. Sampling across a broader species' range would most probably bring about better resolution to differentiate populations based on their metabolite footprints.
In D. ferruginea, the whole population DFPR along with 5 individuals from other 4 populations form a separate cluster ( Figure 6C), while LDA differentiates the population DFZAO from the main cloud formed by the majority of the remaining individuals that belong to 5 other populations. The reason for this differentiation could be the same as for D. grandiflora, namely the limestone soil the plants have been grown on. The most accountable compounds for this differentiation are digoxin, hispidulin, and pcoumaric acid ( Figure 6A) but for the different reasons. Hispidulin in a higher amount is quantified in the population DFPR only, while the other two compounds are the dominant specialized metabolites across the studied D. ferruginea populations, but their quantities vary among these (Gasǐćet al., 2023).
The results point to the fact that the content and ratio of targeted compounds is highly conserved within the analyzed Digitalis species, and is only moderately modified in response to environmental conditions, such as the soil type and altitude above sea level. To increase the resolution of discrimination at the inter-and intrapopulation level, it is possible to include the additional subset of chemomarkers, belonging to the same or other classes of metabolites, which is the course of our further work. Together with morphometrics and molecular markers, chemomarkers can provide a deeper insight into phylogenetic relationships and overall genetic variation among and within Digitalis taxa. Another aspect of this study provides data relevant for the conservation biology studies of the genus Digitalis.
The conservation of individual species, populations, and biodiversity in general is increasingly becoming one of the leading premises of applicative science. The data generated in the present study, which refer to the diversity of bioactive compounds in natural populations of three Digitalis species, thus representing their variation at the phenotypic level, are the starting point for the implementation of conservation strategies. As stated above, digoxin Principal component analysis (PCA) biplot with the first two PCs explaining 83.90% of the total variance (A) among D. lanata accessions. For six D. lanata populations (each labeled by a different symbol) 95% confidence ellipses are presented. Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots. Linear discriminant analysis (LDA) (B) represents the same six populations of D. lanata with 95% confidence ellipses. Heatmap of the scaled quantitative data (C) with the samples arranged according to the hierarchical cluster analysis (Ward's method of cluster agglomeration). Compounds' labels are the same as in Figure 3. For the interpretation of population abbreviations in the figure legend, please refer to Table 1. and other commercially interesting cardenolides are mainly isolated from natural sources, which might severely reduce the natural diversity in wild populations. The logical extension of this research is the study of diversity at the genotype level using various molecular markers, which is also the course of our further work.

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: http://radar.ibiss.bg.ac.rs/ handle/123456789/5448. Principal component analysis (PCA) biplot with the first two PCs explaining 83.84% of the total variance (A) among D. ferruginea accessions. For seven D. ferruginea populations (each labeled by a different symbol) 95% confidence ellipses are presented. Participation of the variables (compounds) in the first two PCs is indicated by the corresponding loading plots. Linear discriminant analysis (LDA) (B) represents the same seven populations of D. ferruginea with 95% confidence ellipses. Heatmap of the scaled quantitative data (C) with the samples arranged according to the hierarchical cluster analysis (Ward's method of cluster agglomeration). Compounds' labels are the same as in Figure 3. For the interpretation of population abbreviations in the figure legend, please refer to Table 1.