Metabolic variation in Cistus monspeliensis L. ecotypes correlated to their plant-fungal interactions

The effect of environmental factors on the chemical composition of plants eventually resulting in plant growth regulation is an age-old issue in plant biology. Nowadays, the acceleration in changes in environmental conditions (e.g. global warming) can act as an incentive to investigate their correlation with metabolic changes. In this study, Cistus monspeliensis plants grown on the island of Sardinia (Italy) were used to explore the geo-graphical-mediated metabolic variation and its repercussion on plant-fungus interactions. Samples of different ecotypes of C. monspeliensis were collected and chemically profiled by 1 H NMR and HPTLC-based metabolomics and the relationship between the variations of biological activity was examined by multivariate data analysis. The ecotypes, collected from different geographical zones and altitudes, exhibited clearly distinguishable chemical profiles, particularly in their terpene and phenolic contents. In particular, multivariate data analysis re- vealed several diterpenes of the labdane and clerodane series among the terpenes and methoxyflavonoids to be responsible for the differentiation. The antifungal activity of the plants was used to explore the correlation between chemical variation and biological activity. Results showed that there was a strong correlation between the metabolic profiles and the antifungal activity, revealing terpenes and methoxylated flavonoids as the main involved metabolites. This demonstrated that environmental factors can influence the chemical variation of plant ecotypes, resulting in the generation of chemotypes that are potentially adapted to their niche conditions including the plant-fungal interactions. to δ 1.10 suggested presence of 8-hydroxylabdan-15-oic acid. confirmed 1 H NMR spectrum of the isolated compound. In the case of the clerodane analogues, the shift of the H-18 chemical shift indicated the presence of α-oxygen and α-olefinic group as in 18-oic acid methyl ester-clerodan-l5-oic acid. This was further confirmed by the presence of a triplet at δ 6.55 ( J = 3.9 Hz, H-3).


Introduction
The effect of environmental factors on the processes of plant life such as growth regulation and adaptation, which are mainly mediated by the alteration of their chemical components, is an age-old issue in plant biology (Scott, 1969). Starting with their germination, plants, as part of a living system, have to contend with challenges posed by very dynamic multifactorial environments, developing strategies that allow them to survive and thrive. Environmental factors are usually classified as biotic and abiotic (Andrew et al., 2010). Abiotic factors are generally physical and/or chemical, often derived from the ambient temperature or water and mineral nutrient availability. Biotic factors relate to microbe, herbivore, or plant to plant interactions. Understanding how for their own survival (Bucharova et al., 2016). These organisms include pollinators, endophytes and even parasitoids (Bucharova et al., 2016). Most of these inter-kingdom interactions require a fast response for survival and these responses have been shaped by long-term coevolutionary processes. It is also clear that the kind and number of these organisms is well differentiated by the geographical location (Abrahamson and Weis, 1997;Leimu et al., 2012).
Considering this, the importance of the geographical influence on the production of plant ecotypes cannot be underestimated. An ecotype is a highly divergent group with phenotypic differences that are too low to be considered a subspecies. This phenomenon can occur in the same geographic region but in different ecological niches. It is possible, also, for similar ecotypes to occur in regions that are distant but share ecological conditions (Mayr, 1999). For example, Bucharova et al. (2016) reported differences in the frequency of seed herbivores associated to diverse flowering phenology in plant ecotypes of Centaurea jacea L. It is thus generally accepted that plants from different origins, identified as ecotypes, largely differ in their interactions with other organisms (Bucharova et al., 2016).
The influence of the geographical location on environmental conditions in which a plant species develops is undeniable and is reflected in various phenotypical features. As a result, the composition of the metabolites believed to be involved in many of the physiological and biological functions in plants is also affected. Considering that ecological interactions are mediated by these chemical signals, the study of the metabolomic profiles of ecotypes from diverse locations provides an ideal opportunity to study the correlation between geographical regions, metabolite variation and biological interactions allowing a better understanding of plant specialized metabolite selection and variation (Salomé-Abarca et al., 2019).
In this regard, the rock-rose (Cistus monspeliensis L., Cistaceae), one of the most common Cistus species, is an interesting model. The chemical variation of several of its metabolite groups, influenced by many environmental factors, has been previously reported. The yield of its essential oils has been reported to vary according to the location and this variation was found to be related to the soil type in which the plants grow (Robles and Garzino, 2000) as well as seasonal variation (Palá-Paúl et al., 2005).
Diterpenes are the major components of rock-rose (Robles and Garzino, 2000). The major components of the leaves of the Jara plant (C. landanifer L.) have been reported to be diterpenes and phenolics (Alías et al., 2012). Interestingly, in this species the specific flavonoid and diterpene content in planta was shown to be correlated with environmental factors such as temperature and water regime, though the relative quantitative contents of flavonoids and diterpenes differed in each situation (Alías et al., 2012).
The rock-rose has been reported to proliferate on the island of Sardinia, Italy. Despite having a very low share of our planet's available land, in general they host a large number of endemic species, thus providing an exclusive level of biodiversity (Tershy et al., 2015;Kier et al., 2009;Weigelt et al., 2013). For example, isolated archipelagos have been recognized as natural laboratories for ecological and evolutionary studies (MacArthur and Wilson, 1967;Grant, 1998). Likewise, the island of Sardinia possesses a great environmental diversity because of varying topographical features on one hand and at the same time is isolated from the influence of exotic species due to its separation from the mainland (Fenu et al., 2014). This specific feature made the island potentially interesting for the study of the influence of environmental factors on the chemical composition of a variety of rock-rose ecotypes growing in different geographical locations and how this affected their biological interactions.
With this in mind, 10 different ecotypes of rock rose, i.e., plants naturally growing at different altitudes on the island of Sardinia, were selected and sampled during the same season and profiled by 1 H NMR and HPTLC-based metabolomics. Additionally, the correlation between their metabolic variation and biological activity, in this case, their antifungal activity, was explored by means of multivariate data analysis. The chemical variation and metabolic interaction among ecotypes were, as expected, reflected in a variation of their antifungal activity as confirmed by biological assays.

Results and discussion
For the investigation of the effects of geographical origin on the plant metabolome, 10 ecotypes of rock-rose (C. monspeliensis L.) were collected from several locations of Sardinia in Italy during the early spring of 2017. The methanol extracts of five specimens per location (n = 50) were firstly analyzed by 1 H NMR. The spectra of the extracts showed a high intensity of signals in upfield regions. Several methyl and methylene protons from terpenoids detected in the δ 0.7 -δ 1.5 range, were assigned to labdane and clerodane like compounds (Fig. 1), i.e., δ 0.81 (s), δ 0.87 (s), δ 0.83 (s) δ 1.10 (s) and one doublet at δ 0.94 (d, J = 6.8 Hz). These resonances were assigned to H-18, H-19, H-20, H-17 and H-16, respectively corresponding to methyl substitutions in labdane analogues (Fig. 1A). Additionally, a doublet at δ 0.77 (d), a singlet at δ 0.80 (s), a doublet at δ 0.97 (d, J = 6.7 Hz), and a singlet at δ 1.11 (s) were respectively assigned to the distinctive methyl H-17, H-20, H-16 and H-19 of clerodane derivatives (Fig. 1B). The detailed structure of a labdanoid was elucidated by additional signals: two double doublets (H-14a and H-14 b) at δ 2.28 (J = 14.8 Hz, 6.4 Hz) and δ 2.11 (J = 14.7 Hz, 7.5 Hz) of the side chain of labdanes as well as the shift of one of the methyl signals to δ 1.10 suggested the presence of 8hydroxylabdan-15-oic acid. This was further confirmed by the 1 H NMR spectrum of the isolated compound. In the case of the clerodane analogues, the shift of the H-18 chemical shift indicated the presence of an α-oxygen and α-olefinic group as in 18-oic acid methyl ester-clerodan-l5-oic acid. This was further confirmed by the presence of a triplet at δ 6.55 (J = 3.9 Hz, H-3).
To visualize the metabolic differentiation of the plants, principal component analysis (PCA) was applied to the 1 H NMR dataset. With eight PCs, the explained variation (R 2 ) reached 0.89. The large variation within the ecotypes, however, interfered in their separation when applied to PC1 and PC2 ( Supplementary Fig. 1A). The ecotype from Su Dominariu (Su) appeared to be distinguished from the other ecotypes, mainly in PC1, though one of its replicates was located far from the cluster. Moreover, despite partial overlapping, the samples from Gennargentu (Ge) and Paringianu (Pa) were distinguished from others in PC2. The samples from Foresta fontanamela (Fo) showed a possible clustering along the PC2.
To supplement the results of 1 H NMR data, an HPTLC analysis targeted on lipophilic terpenoids and phenolic compounds was performed on the same sample set. The samples were analyzed on silica gel plates with a mobile phase consisting of toluene-isopropanol-formic acid (9:1:1, v/v/v). The PCA obtained with the resulting HPTLC data (nine PCs, R 2 = 0.83) did not provide an adequate separation of most of the ecotypes with the exception of the samples from Cardeu (Ca) and Su Dominariu (Su) in PC1 ( Supplementary Fig. 1B).
Although PCA might have the power to extract hidden clustering information in an unsupervised manner, there are two main limitations: the achieved variation is not necessarily related to the clustering, and the available PCs are limited to a maximum of three PCs for plot visualization. Therefore, keeping the unsupervised character of PCA, all the calculated PCs (eight PCs in this case) for the separation were used and a soft independent modelling of class analogy (SIMCA) analysis was chosen as a follow-up analysis. In SIMCA analysis a PCA model is built for each class, and the distances between models (DModX) are measured to determine their degree of similarity. Table 1 show the results of the DModX comparison for 1 H NMR data: the distance of each ecotype in the column was calculated from the averaged distances of the five samples of each ecotype with respect to the set class (ecotype) in the model. When the DModX values of the compared samples from a set class in the model are below 1.5 DCrit (the critical distance to model), they are considered to be clearly classified in their own models. Thus, the DModX values of the same samples in a SIMCA analysis are expected to be below 1.5. Conversely, values above 3 are to be expected for samples that are significantly distinguished from the set class. Thus, the more compared samples with DModX values above 3, the better the set classes were grouped. The compared ecotypes having DModX significantly higher than 3 are indicated in red in Table 1.
According to the SIMCA analysis based on 1 H NMR data, the most distinguishable ecotypes were Fo and Su that showed DModX values above 3 for all the samples. The Se samples, however, were not distinguishable from any of the other ecotypes. This might explain the random distribution of these samples in the PCA score plot. The rest of the ecotypes were partially distinguished from some other groups (Table 1).
The clustering by SIMCA analysis was also performed with the HPTLC data (Table 2). Similarly to the results obtained with the 1 H NMR data, the samples of Fo were clearly discriminated from other ecotypes and Ge samples were also clustered together. However, while the Su sample had been discriminated with 1 H NMR, this did not occur with the HPTLC analysis. Conversely, a greater separation was observed for the samples from Ca and Ma with the HPTLC data. The discrepancy between 1 H NMR and HPTLC could likely be due to the difference in the type of metabolites detected in each case: HPTLC analysis targeted lipophilic terpenoids and phenolics while 1 H NMR detected a wide range of metabolites. It could be thus concluded that the ecotypes are characterized mainly by their terpene and phenolic profiles allowing a greater discrimination based on the HPTLC data. The samples from Barbusi (Ba) as well as Se could not be distinguished from any of the other ecotypes ( Table 2). The SIMCA analysis of both 1 H NMR and HPTLC of the ten analyzed ecotypes, allowed the best differentiation of Fo, Su and Ge while Se and Ba were the least discriminated.
To score the 1 H NMR signals related to the metabolic variation among ecotypes, their standard error was calculated using the buckets of the 1 H NMR spectra. The plot of the average values against their chemical shift revealed the highest variation in the δ 0.6 − δ 2.4 range (Fig. 2), characteristic of terpene analogues. This was confirmed by the PCA loading plot, that showed terpene signals (δ 0.84 -δ 0.92, δ 1.54 − δ 1.68) to be more related to samples from Fo. The loadings from Ge were correlated to terpene olefinic protons affected by the proximity of oxygen (δ 0.76, δ 1.84, δ 1.96, δ 2.20 -δ 2.28, δ 5.68). The Su samples showed signals related to phenolics (δ 7.16 -δ 8.00). Thus, the main differences among ecotypes were due mainly to terpene-like compounds and phenolics.
As SIMCA analysis was unable to yield a clear clustering of all the ecotypes according to their locations, the 1 H NMR and HPTLC data set was further processed using two supervised multivariate data analysis methods, projection to latent structures discriminant analysis (PLS-DA) and an orthogonal projection to latent structures discriminant analysis (OPLS-DA) using geographical locations as categorical variables.
To validate the model, each Q 2 -value was monitored as an indicator of the correlation between the chemical data set and the classes (geographical origins). Only the values above 0.40 were considered as validated (Xiaoying et al., 2011;Cai et al., 2012). Moreover, the analysis of variance testing of cross-validated predictive residuals (CV-ANOVA) was used to assess the reliability of the supervised models. For this test, p-values < 0.05 were considered as significant and p < 0.01 as highly significant. In all, neither PLS-DA nor OPLS-DA models allowed the full separation of all the ecotypes.
The environment in each location of the island of Sardinia includes herbivores, microorganisms, temperature, and rainfall, and many other biotic or abiotic factors that may be unique in some aspects but quite similar in others. Based on this assumption, several OPLS-DA models that showed different groupings of the locations were tested. Interestingly, the model based on regional grouping such as North, Central, and South as categorical classes, was validated. The relative and arbitrary topographical classes were determined according to the closeness of the ecotypes in the sampling area (Fig. 3). When OPLS-DA was applied to the 1 H NMR data, well-clustered groups depending on the locational grouping was achieved after excluding Se samples that had not been distinguished from others at all in the SIMCA analysis (Fig. 4A). The Q 2 -values for each group were 0.61 (Central), 0.55 (South) and 0.5 (North), respectively. Together with CV-ANOVA test (p < 0.01) the model was well validated.
A similar result was obtained with the HPTLC data. Grouping the data according to the regions also revealed a large effect on the metabolome of the plants (Fig. 4B). The validation data, the Q 2 -values of 0.53, 0.45 and 0.80 for the North, South and Central regions respectively with a p < 0.01 for the CV-ANOVA proved the significance of the model (Fig. 4B). Thus, it was possible to conclude that environmental similarities led to these groupings. However, it was not possible to identify which factor(s) contributed to this grouping with the available information. What was clear was that altitude was not relevant in the metabolomic profile of the samples, since there was a wide range of altitudes within each grouped location.
A similar study performed on the essential oils of rock-rose leaves growing in calcareous and siliceous soils in Provence (France), revealed intraspecific chemical variation (Robles and Garzino, 2000). However, the correlation between the observed chemical variation and biotic or abiotic factors remained unclear. In another case, the chemical profile of the essential oils of Cistus parviflorus L. plants collected in nine different locations in the island of Crete showed a large chemical variation among samples with the exception of some mono-and sesquiterpene compounds that were found to be dominant in the samples from all locations . This was not the case for the chemical composition of essential oils from plants of another species, Rosemary (Rosmarinus officinalis) that proved to be influenced by their genetical background rather than their location (Angioni et al., 2004). The distances data are average values (n = 5) calculated from ten rock-rose ecotypes. The values represent the average distance of the ecotypes in the first row (the class marked in italics) from the ecotypes in the first column (the set class marked in bold). When comparing the distances for the same ecotypes, i.e., the distance that figures in the intersection of the corresponding row and column, the values are expected to be below the critical distance to the model (DCrit) value (1.50) and they are thus well classified in their own class (Geographical origin). Those samples with DModX values above 3 are considered to be well-distinguished from the set class in the model. L.F. Salomé-Abarca, et al. Phytochemistry 176 (2020) 112402 The influence of geographical location per se cannot be reduced to one only factor but is more likely the result of a number of interrelated features such as temperature, humidity, rainfalls, soil conditions including its microbiome and altitude among others (Körner, 1999). In this study, one of the factors that seemed most relevant due to its range of variation in the island, was altitude. Thus, the data set was analyzed for the correlation between chemical profiles and altitude. This correlation was examined by OPLS modelling, using the altitude data as the Y variable. A strong correlation was observed with the validation of CV-ANOVA and permutation test as well as a high Q 2value (0.75 for 1 H NMR data and 0.62 for HPTLC data) ( Fig. 4C and D). According to the loading scatter plot and the variable importance for the projection (VIP) plot, one of the signals found to be related to altitude variation was δ 3.12 (dd, J = 7.9 Hz, 1.3 Hz) assigned to one of the H-14 of 8,15-labdanediol. Also, two singlets at δ 2.03 and δ 2.01 were attributed to protons in the methyl ester group of 18-oic acid methyl ester-clerodan-l5-oic acid, and 8-hydroxylabdan-15-oic acid methyl ester, respectively. Furthermore, a singlet in the δ 4.7 − δ 4.6 region was attributed to the H-18 also from 18-oic acid methyl esterclerodan-l5-oic acid and the singlet at δ 9.35 assigned to an aldehyde resonance at C-18 confirmed the compound as an aldehyde derivative (Kalpoutzakis et al., 2002).
Besides the terpenoidal signals detected in the δ 1.9, δ 2.3 and δ 4.0 -δ 5.8 ranges (adjacent to the α-oxygen and olefinic bond), phenolic signals at δ 7.8 − δ 8.4 were also strongly correlated to the altitudinal variation. Phenolics including flavonoids have been depicted as marker metabolites for altitude. The level of phenolics in plants has been reported to be strongly associated to variations in altitude, a phenomenon that was suggested to be explained by their exposure to different levels of UV radiation at each altitude (Alonso-Amelot et al., 2004). For example, environmental factors related to high altitudes resulted in increased flavonol contents in aerial organs of rosebay willow herb (Epilobium angustifolium) (Monschein et al., 2015). On the other hand, this association between altitude and phenolic variation has also been associated with temperature variations as observed for different ecotypes of Arnica montana L. cv. ARBO. In this case, the ratio between the flavonols with an ortho-dihydroxy-configuration in the B-ring (e.g.quercetin) and those with mono-hydroxy groups (e.g. kaempferol) was significantly higher when the plants were grown at lower temperatures (Albert et al., 2009). Moreover, in boxwood (Buxus sempervirens L.) plants grown at different altitudes, the variations in flavonoid contents were found to be related to ontogenic effects rather than UV-B effects (Bernal et al., 2013). In the gum rockrose (Cistus ladanifer L.), the concentration of flavonoids and terpenoids in leaves showed a remarkable seasonal and especially climate-related variation, principally due to temperature fluctuations. While higher temperatures were observed to increase the level of flavonoids, low temperatures increased the amount of diterpenes in leaves (Alías et al., 2012).
It could be assumed that if environmental factors influence the metabolome of plants and these were capable of producing inducible metabolites for their adaptation, there should be a correlation between the metabolome and pathogens, at least to some degree. To prove this hypothesis, two ubiquitous fungal pathogens were chosen to challenge the ecotypes. Fusarium oxysporum was chosen as a soil pathogen and Botrytis cinerea as an aerial pathogen. Strictly speaking, though F. oxysporum is a root pathogen, at some point it can colonize aerial organs through the xylem vessels (Srinivas et al., 2019). A significant antifungal activity of Cistus methanolic extracts was shown only against F. oxysporum. Thus, a PLS and an OPLS model was constructed using the chemical and antifungal activity data obtained in the tests against F. oxysporum. Neither model revealed a clear correlation between the metabolome and the antifungal activity. Yet, considering the result of the biological assay that showed a correlation between the antifungal activities against F. oxysporum and the tested extracts, the Su samples that had shown an exceptionally large intra-set variation, were excluded. The model was now validated with a Q 2 -value that increased from 0.1 to 0.41 in the OPLS analysis (Fig. 5A).
Using the loading and variable importance for the projection (VIP) plots, some of the identifiable loadings related to the antifungal activity against F. oxysporum were deduced to be singlets at δ 0.80 and δ 0.92 assigned to the methyl protons on C-18 and C-16 in a labdane skeleton. This indicated that qualitative and/or quantitative differences in the content of these terpenoids are associated to differences in the antifungal activity among ecotypes. Additional confirmation was provided Gennangertu.

Fig. 2.
Average standard error variation for the buckets of 1 H NMR spectra from ten ecotypes of Cistus monspeliensis L.
L.F. Salomé-Abarca, et al. Phytochemistry 176 (2020) 112402 by loadings of singlets in the δ 2.0 − δ 2.4 range that might be assigned to acetate moieties of 18-oic acid methyl ester-clerodan-l5-oic acid or 8hydroxylabdan-15-oic acid methyl ester, both of which were also associated to differences in antifungal activity. In addition, the singlet at δ 4.59 assigned to the H-18 of 18-oic acid methyl ester-clerodan-15-oic acid was also correlated. Furthermore, the signal at δ 2.11 (dd, J = 14.7 Hz, 7.5 Hz) of H-14 b of 8-hydroxylabdan-15-oic acid also allowed the correlation between variations in the level of this compound and differences in antifungal activity. Besides the signals of these diterpenoids, two characteristic methoxy resonances of flavonoids were found to be correlated with the activity. When the methanolic extract of Ca, which showed higher intensities for 8-hydroxylabdan-15-oic acid signals, was compared with the one from Pa in which this metabolite was not detected, a relatively higher antifungal activity was observed for Ca. This difference was observed as an increase in its inhibition halo of 1 mm. This suggested that the higher content of this diterpene could lead to higher antifungal activity. However, neither of these two ecotypes were among the most active ecotypes ( Supplementary Fig. 2). It could be thus concluded that other metabolites, possibly the methoxy flavonoids mentioned before, might possess a complementary role in the antifungal activity, potentiating the diterpenoids.
Further identification of the metabolites responsible for the antifungal activity was carried out using HPTLC analysis. When HPTLC data of the samples were correlated with the antifungal activity against F. oxysporum, the OPLS model was well validated by Q 2 -value (0.45) and CV-ANOVA (p < 0.01), especially when excluding the Su ecotypes as in the case of the 1 H NMR analysis.
The identification of the correlated metabolites was carried out using a bioautographic assay with HPTLC. This technique has the advantage of allowing compounds in targeted bands to be directly scraped off the plates and isolated. For this, the methanol extracts of the samples of each ecotype were combined to obtain ten representative extracts, separated by HPTLC and tested by bioautography against F. oxysporum. The R f of the bands exhibiting antifungal activity were the same for all ecotypes (Fig. 6). These inhibition zones (R f 0.20 -R f 0.29, R f 0.35 -R f 0.55) were consistent with the data of loading and VIP plots of the OPLS model. The compounds present in the isolated bands were identified by GM-MS as 8,15-labdanediol, 8-hydroxylabdan-15-oic acid, 18-methyl ester-clerodan-15-oic acid, and 8-hydroxylabdan-15-oic acid methyl ester, and confirmed by co-HPTLC with the isolated compounds ( Fig. 7A and B).
Everything considered the presence of both diterpenoids and methoxy flavonoids can be associated to differences in the antifungal activity amongst these ecotypes. Both labdane and clerodane diterpenoids have previously been reported to be among the main bioactive ingredients of the Cistus genus (Papaefthimiou et al., 2014). As for flavonoids, apart from their proven antimicrobial activity, they have been recognized as playing significant roles in plant interactions with microorganisms, including both bacteria and fungi, and their environments (Mierziak et al., 2012). However, the fact that both types of metabolites were present in active ecotypes could support the possibility that the observed antifungal activity is the result of a synergistic interaction between these metabolites. Considering a study that reported the bacterial cell disturbing properties of the flavonol myricetin 3,4′,5′,7-tetramethyl from C. monspeliensis , it could be proposed that this membrane activity could in fact respond to a synergism with diterpenoids.
To explore this possibility of synergistic interaction, the antifungal activity of the isolated compounds against F. oxysporum was first tested individually and then in mixtures. Individually, the most active metabolite was 8,15-labdanediol, with a minimum effective concentration (MEC) value of 2.34 mM after 24 h followed by 18-methyl ester-clerodan-l5-oic acid (MEC = 9.37 mM) and 8-hydroxylabdan-15-oic acid (MEC = 18.75 mM). On the other hand, none of the flavonoids, displayed antifungal activity ( Supplementary Fig. 3). Interestingly, the activity of 8,15-labdanediol disappeared almost completely at all concentrations after 48 h, indicative of a fungistatic type of activity (Banihashemi and Abivardi, 2011). However, the less active 8-hydroxylabdan-15-oic acid still showed full growth inhibition at 150 mM and a strong antifungal effect at 75 mM ( Supplementary Fig. 3). These results showed that rock-rose plants produce highly diverse types of antifungal compounds which vary in strength and duration of their effect, displaying various defense strategies (Bouamamaa et al., 2006). It could also be possible that fungi produce their own metabolites which could counteract the antifungal effect of the tested metabolites, explaining the loss of activity in the microdilution test.
Combinations of the two types of isolated compounds in a 1:1 M ratio were then tested for antifungal activity. The addition of myricetin 3,7,4′,5′-tetramethyl ether to 8-hydroxylabdan-15-oic acid increased its activity, halving its MEC (Supplementary Fig. 3). Hence, this methoxylated flavonoid produced an additive effect on the antifungal effect of this diterpenoid, not having per se a direct antifungal effect as mentioned previously. Another interesting feature of this particular combination was its specificity. The same flavonoid had no effect when combined with the other diterpenoids, at least in the 1:1 M ratio. Thus, the presence of other methoxylated flavonoids and/or other possible molar combinations merit further studies on the rock-rose as an excellent model to explore methoxylated flavonoid-diterpenoid interactions not limited to additive effects but also potentially synergistic effects.

Conclusions
The metabolome of rock-rose (C. monspeliensis L.) was influenced greatly by topographical and altitudinal variations in the island of Sardinia. Of the metabolites, diterpenoids and methoxy flavonoids were the compounds correlated to the ecotypes' metabolic differentiation. One of the premises of this study was that metabolic variation in plants have developed as a response to environmental factors should inevitably be reflected on their activity against pathogens, given that they are part of the environment with which the plants interact. This in fact proved to be true for the systems analyzed in this study since the observed chemical variation among the tested ecotypes was well correlated with the activity against the fungus F. oxysporum. Based on the activity test of several isolated diterpenoids and flavonoids found to be involved in the antifungal activity it was possible to reach some conclusions on the adaptation of rock-rose to its environment. Being confined to the island of Sardinia, the diverse Cistus ecotypes are a result of an adaptation to geographical variation, producing diverse types of metabolites with different characteristics regarding their antifungal L.F. Salomé-Abarca, et al. Phytochemistry 176 (2020) 112402 effects exhibiting long-lasting, instantaneous, direct, supplementary or indirect effects according to the environmental challenges of each location.

Plant collection
The

1 H NMR analysis
Plant material (30 mg) was extracted in 1 mL of CH 3 OH-d 4 containing 3.93 mM hexamethyldisiloxane (HMDSO) as an internal standard followed by 20 min of ultra-sonication. The extracts were centrifuged at 13 000 rpm for 10 min, and 300 μL aliquots of the solutions were transferred to 3 mm-NMR tubes for 1 H NMR analysis. The 1 H NMR analysis was carried out in an AV-600 MHz NMR spectrometer (Bruker, Karlsruhe, Germany), operating at proton NMR frequency of 600.13 MHz. For internal locking CH 3 OH-d 4 was used. All 1 H NMR analyses consisted of 64 scans requiring 10 min and 26s as acquisition time using the parameters: 0.16 Hz/point, pulse width (PW) = 30°( 11.3 μs), and relaxation time of 1.5 s. A pre-saturation sequence was used to suppress the residual water signal using low power selective irradiation at H 2 O frequency during the recycle delay. The FIDs were Fourier transformed with exponential line broadening of 0.3 Hz. The resulting spectra were manually phased and baseline corrected, and calibrated to HMDSO at δ 0.06 using TOPSPIN V. 3.0 (Bruker).

Compound isolation and identification
Fifty g of dried plant material of each of the ecotypes Mandas (Ma) and Seui (Se) were extracted with 500 mL of MeOH-H 2 O (1:1, v/v) with ultra-sonication for 1 h. The extracts were filtered and taken to total dryness with a rotary evaporator yielding residues of 33% (w/w) and 31% (w/w) for Ma and Se respectively. These residues were redissolved in 10% aqueous methanol and partitioned with n-hexane.

High performance thin layer chromatography (HPTLC)
Silica gel HPTLC plates (20 × 10 cm, F 254 ) were purchased from Merck (Darmstadt, Germany). Samples were processed using a CAMAG-HPTLC system equipped with an automatic TLC sampler (ATS4), a derivatizer (version 1.0 AT), a TLC plate heater (version III), a TLC scanner (version 3) and a TLC visualizer (version 2) (CAMAG, Muttenz, Switzerland). For sample extraction, 30 mg of plant material were extracted with 1 mL of methanol and centrifuged at 13 000 rpm for 10 min. The supernatant was used for the analysis. Ten μL of each sample were spotted as 6 mm bands. Fifteen samples were spotted on each plate leaving 10 mm distance from the bottom and 20 mm from the left and right edges of the plate. The distance between bands was 11.4 mm. The distances were set to be calculated from the center of the band. A pool of all samples was applied at the right end of each plate as a quality control sample (QC). The mobile phase for HPTLC analysis was a mixture of toluene-isopropanol-formic acid (9:1:1, v/v/v). Saturation time of the chamber was 20 min and the solvent migration distance was 80 mm from the application point. The developed HPTLC plates were sprayed with 2 mL of anisaldehyde-sulfuric acid solution using an automatic derivatizer and placed on a TLC plate heater at 100°C for 3 min. Images of the plates were recorded using a TLC visualizer at 254 and 366 nm before derivatization, and 366 nm and white light after derivatization.
When the scanner was used for band spectrum acquisition, it was set to automatic band detection at 254 nm and 366 nm, then the chromatographic bands were transformed into densitograms. After band assignment, the scanner was used to acquire a UV/Vis absorption spectrum between 200 nm and 500 nm. The spectrum acquisition speed was 20 nm/s and the band slit was 4 × 0.1 mm. The whole instrument was controlled using the software VisionCats (V. 2.5, CAMAG).

GC-MS analysis
The bands scraped from the HPTLC plates were ultra-sonicated with 1 mL of acetone for 3 min, centrifuged at 13 000 rpm for 10 min to separate the supernatants which were removed and taken to dryness with a Speedvac. The resulting residues and the pure isolated compounds were then derivatized by adding 100 μL of pyridine (Sigma-Aldrich, Zwijndrecht, The Netherlands) and 100 μL of BSTFA + TMCS (99:1) (Sigma-Aldrich, Zwijndrecht, The Netherlands) and heating at 80°C for 50 min. The GC-MS analysis was performed using a 7890 A gas chromatograph equipped with a 7693 automatic sampler and a 5975C single-quadrupole mass detector (Agilent, Folsom, CA, USA). Samples were separated on a DB-5 GC column (30 m × 0.25 mm, 0.25 μm film, J&W Science, Folsom, CA, USA) and eluted with Helium (99.9% purity) as a carrier gas at a flow rate of 1 mL/min. The oven temperature was programmed as follows: initial temperature was held at 60°C for 1 min, then increased at 7°C/min to 290°C and held for 5 min. The injector was set at 280°C and 1 μL of sample was injected with a 5:1 split and split flow of 5 mL/min. The interface temperature was 280°C, and the ion source and quadrupole temperature of the mass detector were 230°C and 150°C, respectively. The ionization energy in EI mode was 70 eV. Compounds were identified by comparison of their retention times and ion spectra with those of the isolated compounds. Data was processed using Mass Hunter (B.07, Agilent), AMDIS (V. 2.62, Agilent, Folsom, CA, USA) and MS search (V. 2.0, Agilent).

Data processing and multivariate data analysis
The NMR spectra were bucketed using AMIX 3.9.12 (Bruker BioSpin, Rheinstetten, Germany). Bucket data was obtained by spectral integration at 0.04 ppm intervals. Peak intensity of individual peaks was scaled to the total intensity recorded from δ 0.20 to δ 10.00. Because of the residual signals of D 2 O and CH 3 OH-d 4 , regions of δ 4.75 -δ 4.9 and δ 3.28 -δ 3.34 were excluded from the analysis, respectively.
The software rTLC (version 1.0) was used for data extraction and alignment of HPTLC images, according to Fichou and colleagues (2016). The same dimensions were used for image analysis and sample application on the silica plates. For data binning, 128 units (pixel width) were employed. A parametric time-warping algorithm was applied to the R f data for band alignment. The data from the grey channel was used for multivariate data analysis (MVDA). After data processing, all signals were normalized to the corresponding pooled QC sample signals of each plate.
To observe the variation of the secondary metabolites across all ecotypes, the standard error of each bucket was calculated and plotted against its corresponding chemical shift using Microsoft Excel 2013. Multivariate data analysis was done with SIMCA P (v.15.0.2, Umetrics, Umeå, Sweden). Both 1 H NMR and HPTLC data were submitted to principal component analysis (PCA), soft independent modelling of class analogy (SIMCA) and orthogonal partial least squares (OPLS). The data was scaled by unit variance (UV) scaling method after mean-centering for both 1 H NMR and HPTLC. In the SIMCA analyses, a PCA model was built for each class and the distance between models (DModX) was measured to determine their similarity or dissimilarity. The distance of each ecotype was calculated from the averaged distances of the five samples of each ecotype to the set class (ecotype) in the model. When the values were below the critical limit for the distance to model (DCrit), that is, 1.5, the samples were interpreted to be non-discriminated from the set class. Thus, the distances for samples belonging to the same cluster, geographical origin in this case, would all be expected to be below the DCrit. For comparison purposes, the samples from other classes with DModX values above the DCrit value but lower than 3 were considered to be similar (non-discriminated) to the set class in the model. The DModX values of all ecotypes based on 1 H NMR and HPTLC data were tabulated and presented as averaged L.F. Salomé-Abarca, et al. Phytochemistry 176 (2020) 112402 DModX values (n = 5) ± standard deviation. Models for OPLS-DA and OPLS analyses were validated by using the Q 2 -value as an indicator of the correlation between the chemical data set and the classes (geographical origins) according to a cross-validation test. Negative Q 2 values indicated no correlation between the X matrix and the set classes, proving a total lack of predictability of the model. Positive Q 2 values indicate a degree of correlation, but in this case only values of Q 2 ≥ 0.40 were considered to be validated (Xiaoying et al., 2011;Cai et al., 2012). The Analysis of variance testing of Cross-Validated predictive residuals (CV-ANOVA) was used to gauge the reliability of the supervised models. For this test, p-values < 0.05 were considered to be significant, and p < 0.01 as highly significant while values between 0.05 and 0.10, that is, those indicating a 5%-10% of probability of not getting the expected results, were considered to be indicators of data trends.

Medium composition and preparation
The complete medium consisted of 10 g of D-glucose, 5.96 g of NaNO 3 , 1.49 g of KCl, 0.5 g of MgSO 4 .7H 2 O, 1 g of casaminoacids, and 5 g of yeast extract dissolved in a 2% aqueous solution of agar. The medium was sterilized at 121°C for 20 min. The 2X minimal medium was an aqueous solution containing 10 mL of ASPA + N ((50X) prepared by dissolving 29.75 g of NaNO 3 , 2.61 g of KCl, and 7.48 g of KH 2 PO 4 in 60 mL MilliQ water, adjusting to pH to 5.5 with KOH and taking to a final volume of 100 mL with MilliQ water), 10 mL of 50% Dglucose solution, 1 mL of 1 M MgSO 4 , 0.5 mL of 1000x Vishniac, and 5 mL of a 0.3% yeast extract, taken to 250 mL with MilliQ water. The medium was sterilized by filtration.

Spore production
The ten methanolic extracts of the ecotypes were tested against Fusarium oxysporum and Botrytis cinerea. To produce F. oxysporum and B. cinerea spores, complete medium plates were inoculated with fungal solutions and incubated at 28°C and 26°C respectively for two weeks in the dark. After this, the plates were filled with 25 mL of sterile physiological-salt solution (Ps) and rubbed with a sterile cotton swab to transfer the mycelium and spores to the Ps. This liquid containing the fungal structures (mycelium and spores) was removed with a sterile pipet and filtered through two layers of sterile miracloth, and transferred to sterile 50 mL conical tubes. The volume was adjusted to 30 mL with sterile water, vortexed for 15 s and then centrifuged at 4000 rpm for 10 min (3x). The supernatant was discarded and the pellet was resuspended in 30 mL of sterile Ps. The spore concentration per milliliter was quantified using a cell counter apparatus (TC20™, Bio Rad, The Netherlands).

Sample preparation
To test for antifungal activity, methanolic extracts of each plant were prepared by sonicating 70 mg of plant material with 1 mL of methanol for 20 min. The mixtures were then centrifuged at 13 000 rpm for 10 min, and the supernatant was removed for the bioassay. Before sample application, silica gel TLC plates (10 × 10 cm, F 254 ) were flushed with ethanol. The excess of solvent was dried with a cold air stream. The extracts were then applied on the plates in five rows (five samples per row) using an automatic TLC sampler version 4 (CAMAG). The distances for spot application were 20 mm on the X axis and 8, 28, 48, 68 and 88 mm in the Y axis. The application volume was 30 μL, and the sample application mode was spray area on square form with 3 × 3 mm as area dimension using a 100 μL Hamilton syringe.

Plate assembly
The bioassay was performed using 12 × 12 cm square petri dishes. The spore solutions were mixed with complete media (50°C) to reach a final concentration of 2.5 × 10 5 spore/mL. Fifty mL of the medium were poured into an empty Petri dish. After medium solidification, the TLC plate was faced up on top of the lower layer and 50 mL of the medium containing the spores were poured over the TLC plate. This resulted in an approximately 80 mm-thick agar layer. After medium solidification, the plates were incubated at 28°C for F. oxisporum and 26°C for B. cinerea in the dark. The inhibition halos were measured after 39 h of incubation.

Minimal effective concentration and potentiation test
The minimal effective concentration (MEC) values of individual compounds and their potentiation were determined by the microdilution and checkerboard assay respectively, based on the Clinical and Laboratory Standard Institute guidelines (NCCLS, 2003) and on previously reported methods (Abreu et al. 2015(Abreu et al. , 2016Bonapace et al., 2000) with some modifications. Briefly, 100 μL Fusarium oxysporum spores were inoculated into 2X minimal medium contained in 96 well plates. The final spore concentration in the well was 2.5 × 10 5 spores/ mL. To reach the germination state, spores were incubated for 6 h at 28°C in the dark. Once the spores germinated, 100 μL of the compounds were added in a mixture of medium and DMSO. The final concentration of DMSO in the well was 1% (v/v) in a final volume of 200 μL. The plates were sealed with a parafilm layer and then incubated in dark conditions at 28°C for 24 h and 48 h. To avoid the interference in the OD readings posed by the precipitation of compounds over time (instantaneously or gradual compound precipitation), microscopic observation at 4X and 40X of the wells was used to check the fungal growth. The MEC was defined as the lowest drug concentration that visually decreased the total fungal growth with respect to the negative control after 24 h and 48 h. The treatments consisted in single compound dilutions in the range of 150-15 mM 8,15-labdanediol, 8-hydroxylabdan-15-oic acid, 18-methyl ester-clerodan-l5-oic acid and myricetin 3,7,4′,5′-tetramethyl ether. For the checkerboard assay, combinations of two diterpenes and of each diterpene with the methoxylated flavonoid in a 1:1 M ratio were tested in the same range of concentrations. The negative controls consisted of spores of F. oxysporum in the medium and in the medium with 1% of DMSO (v/v). The positive controls consisted of nystatin in a range of 250 μg/mL to 0.48 μg/mL diluted in sterile water and 0.1% DMSO (v/v). The treatments were randomly distributed in the 96 well plates and the wells of the edges of the plates were not used for compound testing, but instead were filled with the medium and fungal spores. All of the MEC determinations were done by triplicate and checked with a compound microscope.
Metabolite (A) + metabolite (B) interactions were classified by the fractional inhibitory concentration (FIC) index: FICI = FIC(A) + FIC(B), where FIC(A) is the quotient between the MEC of the metabolite A combined with metabolite B and the MEC of metabolite A alone. FIC(B) is the quotient of the MEC of drug B combined with metabolite A and the MEC of metabolite B alone. A FICI value of ≤0.5 was interpreted as synergy, FICI >4 as antagonism and FICI >0.5 and <4 as indifferent (Abreu et al., 2016;Bonapace et al., 2000). For the blends in which one of the metabolites had no antifungal activity, the chemical interactions were defined as synergism when a MEC reduction was observed with 4-fold dilutions, as an additive interaction with a MEC reduction of 2 or 3-folds, and as an antagonist interaction with a MEC increase at 4-fold dilutions (Abreu et al., 2015). Three replicates were performed for potentiation experiments.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.