Plant morphological variability at the distribution edges: the case of Dryas octopetala (Rosaceae) in northern Greece

Plant populations at their distribution limits may often deviate morphologically from those at the centre of their range (core populations). A similar pattern was observed in Dryas octopetala (Rosaceae), a circumpolar, arcticalpine species distributed in northern and central Europe, reaching south to northern Greece. The three southernmost populations of Europe, located in Greece, were sampled and specific leaf morphometric traits were measured and analysed using canonical variates analysis (CVA) and hierarchical cluster analysis. CVA revealed that the three samples of D. octopetala did not form one group but were significantly discriminated. Despite the rather similar climatic conditions of Mts Falakron and Orvilos, the sample of Mt Falakron was grouped with that of Mt Tzena, based on leaf morphology. These findings contradict findings from other studies on the ways with which plants react and adapt on areas of harsh climatic conditions. Micro-site conditions or a potentially different post-glacial origin could possibly explain the pattern observed in this study.


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 268 Varsamis & al.: Morphological variability of Dryas octopetala in northern Greece 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 lim-its 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).
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  Downloaded From: https://bioone.org/journals/Willdenowia on 26 Feb 2021 Terms of Use: https://bioone.org/terms-of-use Willdenowia 50 -2020 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).

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 2 (de Heredia & al. 2009), where π: Archimedes' constant (π = 3.14159) LA: leaf area and LP: leaf perimeter 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 per-formed, 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. 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 negativelycorrelated 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 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 oc topetala. 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 individu-als 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 nonsignificant 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,

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 morphologi-  (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 octo petala 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 Dry as 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 Pen insula (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.