Introduction

Species, as biological organisms, are morphologically variable to a greater or lesser extent. Although there are plant species that are characterized by a restricted morphological variability, especially those that are range restricted or self-pollinated, others can be greatly differentiated throughout their distributional range (e.g. Christensen 1992; Marcysiak 2012; Forsman & Wennersten 2015; Tsiftsis 2016; Zarei & al. 2019).

It is well known that, especially in species with wide distributions, populations occurring at the edges of their range may often deviate morphologically from the more central populations and they also present great morphological variability (Antonovics 1976; Garcia-Ramos & Kirkpatrick 1997; Jonas & Geber 1999). However, specific morphological traits may also be influenced by several environmental factors, e.g. climate, soil (Mal & Lovett-Doust 2005; Peppe & al. 2011; Blinova 2012). Moreover, post-glacial migration and natural selection of specific morphotypes could differentiate fragmented and isolated populations (Hatziskakis & al. 2011; Soubani & al. 2015).

Plants that exhibit large intraspecific variation (e.g. in physiology, morphology, phenology) constitute ideal cases for studying local or regional adaptations (Gratani 2014). Until now, several studies have been published exploring the range of morphological variability of species or referring to plant discrimination. Such studies use a variety of plant traits, mostly depending on the taxonomic group (family or genus) to which the species concerned belongs. In some studies, flower morphology has been used to explore the morphological variability of different populations (e.g. Medrano & al. 2006; da Cruz & al. 2011; Tsiftsis 2016). Leaf morphology has also been widely used as a taxonomic tool by botanists, mostly for separating taxa at the specific or even subspecific level (e.g. Viscosi & al. 2009; de Oliveira & al. 2018).

Leaf morphology is widely considered as a group of trait variables from which biogeographical information and species evolution can be deduced (Klingenberg 2010; Dkhar & Pareek 2014; Marcysiak 2014; Tian & al. 2016; Lai & al. 2018). This is because leaves are essential parts of the adaptation of a species to a given ecological niche, and such information can be expressed through a number of leaf morphological traits. A previous study demonstrated that both leaf size and leaf shape of Dryas octopetala L. (Rosaceae) follow a biogeographical pattern (Marcysiak 2014). Specifically, it was found that D. octopetala populations are clustered into several distinct groups, in which the leaf morphological traits (e.g. leaf area, number of teeth, maximal width of the apical tooth) are mainly influenced by climate, throughout Europe. Based on these findings, Marcysiak (2014) considered D. octopetala as an ideal species for exploring geographical patterns. Moreover, the study of widespread species at their southernmost distribution limits is interesting because it can provide valuable information about their ecology (Schuler 2004).

Dryas octopetala (Fig. 1) is a circumpolar, arctic-alpine plant species growing in areas where climatic conditions are harsh. In the north of its distributional range, it can be found even in low-altitude areas, whereas at its southern distribution limits it exclusively occurs high on mountains (Tutin & al. 1968; Strid 1986). It is a woody, long-lived dwarf shrub (life expectancy over 500 years; de Witte & al. 2012) and constitutes the ground vegetation cover in the arctic and alpine tundra ecosystems. The populations of the species currently extend from Eurasian and American arctic tundra to the temperate mountains of Europe, including the Balkan Peninsula (Hultén & Fries, 1986).

Fig. 1.

Dryas octopetala in flower. – Greece, Eastern Macedonia, near summit of Mt Falakron, 26 June 2006, photograph by S. Tsiftsis.

In her study, Marcysiak (2014) sampled material of Dryas octopetala from numerous populations, almost throughout the European range of the species, although not from northern Greece, its southernmost distribution limit in Europe. In Greece, D. octopetala is considered a remnant species, which probably migrated there during the last glacial age, now forming three distinct and isolated populations in the northern part of the country (Schuler & Tsiripidis 2009). The fact that the Greek populations of D. octopetala have been fragmented and isolated for a long time (maybe for several centuries or even millennia) makes them an ideal case to study their morphological differentiation or divergence. Based on this background, the aim of the present study is to examine whether there are leaf morphological differences between individuals from the three Greek populations, and if any differentiation could be explained by the biogeographical patterns identified by Marcysiak (2014).

Fig. 2.

Map of the studied populations of Dryas octopetala.

Material and methods

We collected leaf material from the three known Dryas octopetala populations in Greece, namely Mt Tzena (41°09′N, 22°14′E), Mt Orvilos (41°22′N, 23°37′E) and Mt Falakron (41°17′N, 24°05′E) (Fig. 2). Collection of leaves from individuals of all populations took place in July of 2018. We sampled 35, 22 and 37 individuals, respectively, from each population. A total of 10 leaves was collected from each individual and stored in portable cooler within plastic bags. Leaves were taken from the same height and from the outer part or the plant canopy. The fresh leaf material was scanned and images were processed in black and white mode before morphometric analysis, using GIMP v2.9 software. Following Marcysiak (2014), we measured the following leaf morphometric variables (leaf morphological traits): leaf area (LA), leaf perimeter (LP), leaf circularity (LC), leaf length (LL), leaf width (LW), number of lobes (NL), maximum width of apical tooth (AW), length from midrib to leaf periphery (BL) and leaf base width (BW) (Fig. 3). Leaf circularity, or else called shape factor, was calculated according to the following formula:

Leaf circularity = (4 × π × LA) / LP² (de Heredia & al. 2009), where

Leaf circularity values close to 1 (a circle) indicate more circular leaves compared to values close to zero, which indicate non-circular (more elongated) leaves. All measurements were performed with the freely available software ImageJ (Rueden & al. 2017) and are reported in millimetres (mm). Leaf variables were expressed on an individual plant basis by calculating the average value from 10 leaves per plant.

To study the multivariate patterns of size variation of the measured Dryas octopetala individuals, a Canonical Variates Analysis (CVA) (Mardia & al. 1979) was performed, having the sampled population as the categorical variable. CVA is a widely used technique for assessing and displaying variation among pre-defined groups relative to the variation within the groups. The method successively extracts axes (Canonical Variates, CVs), which are a linear combination of the original variables and orthogonal to all others, having however the greatest ratio of among-group to within-group variance. Therefore, plots of the first few CVs are optimal displays of differences among groups and are scaled relative to the pooled estimate of within-group variation (Webster & Sheets 2010). A cross-validation using the leave-one-out procedure (jackknife test) was performed to estimate the expected actual error rates in classifying D. octopetala individuals in the respective groups.

Fig. 3.

Way of measuring leaf morphological traits of Dryas octopetala; LL: leaf length; LW: leaf width; AW: maximum width of apical tooth; BL: length from midrib to leaf periphery; BW: leaf base width.

Fig. 4.

Representation of the scores on the first two canonical axes of the CVA using nine morphological leaf characters in the three populations of Dryas octopetala.

The morphometric relationships among the three Dryas octopetala populations were further illustrated by a hierarchical cluster analysis using the unweighted pair group method with an arithmetic mean (UPGMA). Cluster analysis was based on the Mahalanobis distances between the means of the samples, as these were calculated by the canonical variate analysis. To test for statistical differentiation between the samples of the three populations, the leaf morphometric data were firstly checked for normality using the Shapiro-Wilk normality test, because this test outperforms others (e.g. Kolmogorov-Smirnov), especially when the sample size is small. After checking for normality, and because not all variables were normally distributed, the differences between the samples of the three populations were tested using Wilk's lambda statistic, calculated using a nonparametric comparison of multivariate samples. The Mann-Whitney U test was used to test for statistical differences among samples of the three populations in each of the nine leaf morphometric variables. Finally, Pearson correlation coefficients between the nine leaf morphological traits and the bioclimatic/climatic variables (obtained from the WorldClim database; Fick & Hijmans 2017) of the three populations were calculated.

All analyses were performed in R version 3.5.2 (R Foundation for Statistical Computing) using the Morpho (Schlager & al. 2019) and npmv (Burchett & Ellis 2017) packages.

Results

The result of the canonical variates analysis is graphically presented in Fig. 4. The first CV axis accumulated 56.2% of the variance and was highly correlated with leaf circularity (LC) and negatively correlated with maximum width of the apical tooth (AW), whereas the second axis accounted for the remaining 43.8% of the total variance and was found to be highly – albeit negatively – correlated with AW and positively with LC (Table 1). Leaf area (LA) was found to be the least significant variable, followed by the number of lobes (NL).

Based on the nine leaf morphometric variables, CVA discriminated the three sampled populations, because they overlap only to a rather small extent. Furthermore, Wilk's lambda statistic demonstrated that the samples of the three populations were significantly discriminated (p < 0.001), whereas the correct classifications for the whole dataset was 70.21%.

The number of correctly classified individuals in the three studied populations is shown in Table 2. The overall percentage of correctly classified individuals was rather high, ranging from 63.63% to 77.14%. The lowest percentage of correctly classified individuals was calculated for the population of Mt Orvilos, whereas the highest was for the population of Mt Tzena. Moreover, it is worth mentioning that eight out of the 12 not correctly classified individuals of Mt Falakron were assigned to the population of Mt Tzena, demonstrating their rather close morphological relationship.

Table 1.

Canonical discriminant coefficients showing the contribution of the nine leaf morphological traits to the two canonical variates (CV) analysis axes.

Fig. 5.

Results of hierarchical cluster analysis (UPGMA) of the three Dryas octopetala populations based on Mahalanobis distances.

Table 2.

Cross-validated (= jackknifed) confusion matrix of the canonical variates analysis (CVA) performed on leaf morphometric data of three Dryas octopetala populations (Mts Falakron, Orvilos and Tzena). The total for each row represents individuals sampled in each population, whereas the numbers in the rows represent individuals classified, on the basis of their morphological affinities, in each population. The numbers in the columns represent individuals classified in each population.

Table 3.

Results of the Mann-Whitney U test comparison among all pairs of mountains. The mountain name in a cell indicates the mountain having a statistically higher value in the specific leaf morphological trait. For trait codes see Table 1; *** p < 0.001; ** p < 0.01; * p < 0.05; ns = non significant.

The UPGMA cluster based on Mahalanobis distances among the pairs of the three sampled populations is shown in Fig. 5. As can be seen, the population of Mt Tzena was grouped together with that of Mt Falakron, whereas that of Mt Orvilos was separate. The greater morphological affinity of the Dryas octopetala individuals from Mts Tzena and Falakron, compared to those from Mt Orvilos, confirmed the greater overlap of these two populations on the CVA graph and the results presented on Table 2 for the correctly classified individuals.

Finally, the differences in the measured leaf traits of the samples of the studied populations are shown in Fig. 6, whereas the results of the Mann-Whitney U test for their statistical differences are presented in Table 3. Specifically, individuals from Mt Falakron are characterized by larger leaves (in terms of area, perimeter, width, length and number of lobes) compared to those from the other two populations, whereas those from Mt Tzena are the smallest.

The pair-wise differences of the populations of Mts Falakron and Tzena were statistically significant for all the leaf traits, except for AW and BW. Contrary to the comparison of these two mountains, in which most differences were statistically significant, the comparison between the other pairs of mountains revealed more non-significant pair-wise differences. Specifically, the comparison between Mts Falakron and Orvilos revealed non-significant differences for LA, LC, LL, BL and BW, whereas the comparison between Mts Orvilos and Tzena revealed non-significant differences for LL, LW, NL, BL and BW.

The Pearson correlation coefficients between the morphological traits and the bioclimatic/climatic variables are shown in Table 4. As can be seen, the effect of both bioclimatic and climatic factors on the morphological traits is in general significant, either positively or negatively. Annual precipitation (bio12) and solar radiation during the summer period (May–September) were factors that significantly affected only one trait, the maximum width of the apical tooth (AW). The number of lobes (NL) and leaf base width (BW) were the least affected morphological traits, as they were affected only by specific bioclimatic and climatic factors (i.e. bio7, bio13, bio16, bio19, specific monthly values of solar radiation and August precipitation), but even in these cases the calculated correlation coefficients were low (< |0.25|). Moreover, harsh climatic conditions, as expressed by the minimum monthly temperatures, negatively affect leaf size (leaf area, perimeter, length and width) and maximum width of the apical tooth, whereas only leaf circularity was positively affected.

Table 4.

Pearson correlation coefficients between bioclimatic (bio) and climatic (srad, prec, tmin) variables with leaf morphological trait variables. For trait codes see Table 1; for bioclimatic/climatic variable codes see Appendix; ** p<0.01; * p<0.05. Values for bioclimatic/climatic variables were downloaded from the WorldClim database (Fick & Hijmans 2017).

Fig. 6.

Boxplots presenting the variation in the nine leaf morphometric variables of the three populations of Dryas octopetala occurring in Greece.

Discussion

The results of the present study confirmed that morphological traits may vary and differ when studying populations at the edges of the distribution of a species (Jonas & Geber 1999; Tsiftsis 2016). Despite the fact that the samples of the three studied populations overlapped to some extent in the CVA graph, they proved to be statistically discriminated. This suggested that either ecological factors (e.g. climatic, soil/geological) or natural selection and a potentially different post-glacial history might have shaped this morphological differentiation.

Dryas octopetala exclusively occurs in sites with limestone, the Greek populations being no exception (Strid 1986; Schuler & Tsiripidis 2009; Marcysiak 2014). On-site examination showed that D. octopetala occurs in Greece above 1900 m in rocky micro-sites of northern or northwestern exposure (see Appendix), where competition with other plant groups (e.g. grasses) is greatly reduced or does not exist, on all three mountains.

Apart from soil/geological conditions and stress due to competition with other plants, climatic conditions have also been found to affect plant morphological traits, including leaf size (Sultan 2000; Peppe & al. 2011). Despite the general worldwide trend of decreasing leaf size with latitude and elevation (Wright & al. 2017), studies that contradict these patterns are not rare in the literature. Specifically, Royer & al. (2008) found that Acer rubrum L. tends to have larger leaves toward colder climates, whereas Quercus kelloggii Newb. has smaller leaves. Based on the results of Marcysiak (2014), Dryas octopetala does not follow the general worldwide trend in its European range. Specifically, Marcysiak (2014) found that leaf size (expressed by LA, LL and BW) was positively correlated with the minimum temperatures during the first part of the growing season (April–June), demonstrating that it increased towards southern regions. Moreover, she found that AW was positively correlated with rainfall during August and September, as well as with the mean temperature of the wettest quarter.

In the present study, we found that leaf size (in terms of LA, LL, LW and LP) is negatively correlated with the minimum temperatures, whereas it is positively correlated with August precipitation but negatively with September precipitation. Moreover, we found that individuals from Mt Falakron had larger leaves compared to individuals from the other two mountain ranges. This finding contradicts the pattern described by Marcysiak (2014), because Mt Tzena has milder climatic conditions (based on data downloaded from the WorldClim database [Fick & Hijmans 2017]; see Appendix), whereas Mt Orvilos has harsher conditions compared to the others. Moreover, under harsh climate conditions, plant leaves tend to have larger lobes and more teeth, and this constitutes an adaptation for increased carbon uptake at the beginning of the growing season (Royer & Wilf 2006; Peppe & al. 2011). However, neither Marcysiak (2014) nor we found such a correlation. Specifically, Marcysiak (2014) did not find any significant correlation among the factors “number of teeth” and “monthly minimum temperatures”, whereas, in our study, the maximum number of lobes (equivalent to the number of teeth sensu Marcysiak 2014) was recorded in the individuals from Mt Falakron.

Besides the differences between our findings and those of other studies, another interesting pattern was the clustering of the Mt Falakron population with the Mt Tzena population, although the former is located closer to the Mt Orvilos population. Mt Falakron and Mt Orvilos are located in eastern Macedonia, forming a continuation of the Bulgarian Pirin Planina (Strid 1986). In that sense, one would expect that the morphological traits of these two populations would be similar and therefore clustered together. However, the observed pattern might indicate an adaptive morphological differentiation, probably due to potentially different micro-ecological conditions that were not analysed in this study. Thus, we could hypothesize that differentiation in leaf traits could be attributed either to different microclimatic conditions, which could not be expressed through the climatic and bioclimatic variables of the WorldClim database (Fick & Hijmans 2017), or to divergent selection and genetic drift caused by the glacial history of the populations. This conclusion was also proposed by Marcysiak (2014) for isolated Dryas octopetala populations, where leaf morphology differentiation was not consistent with the biogeographical molecular analyses conducted by Skrede & al. (2006).

Finally, apart from the effect of microclimate, the studied populations could also represent different Dryas octopetala morphotypes related to a potentially different post-glacial origin. During the last glacial period, in the late Quaternary era, D. octopetala, as many other alpines, migrated toward southern Europe including the Balkan Peninsula (Birks 2008). This upward-downward movement has taken place more than once and caused population fragmentation and isolation after the end of the last glacial episode (Allen & al. 2015). Hence, genetic analyses should be conducted to enable the attribution of the observed differences in morphological traits among the three D. octopetala populations in Greece either to microclimate or to biogeographical patterns of post-glacial origin.