Population Variability of Almond-Leaved Willow (Salix triandra L.) Based on the Leaf Morphometry: Isolation by Distance and Environment Explain Phenotypic Diversity

Abstract

Almond-leaved willow (Salix triandra L., Salicaceae) is a dioecious shrub, rarely a small tree that grows under various environmental conditions. We examined the population structure of 12 populations of almond-leaved willow using nine leaf morphological traits and specific leaf area. Populations were selected from a range of habitats, from continental to the sub-Mediterranean zone, to examine the influence of environmental conditions (climate and altitude) and geographic distance on leaf variability. Significant differences were confirmed among all populations for all traits, with significant correlations between geographic location of populations and morphological traits, and between environmental conditions and morphological traits. Large-leaved populations were found in continental and sub-Mediterranean climates, while small-leaved populations were found in higher elevations and smaller karstic rivers. In addition, populations from floodplains showed greater variability than populations from the karstic habitats, indicating a positive influence of lowland habitats and possible underlying differences in gene pool size. In conclusion, we found that environmental conditions and geographical distances in addition to genetic drift, are the main influences on the variability in almond-leaved willow, with the species showing a high level of plasticity and adaptation to local environmental conditions.

1. Introduction

Leaves are essential organs for photosynthesis and carbon assimilation [1], making their morphological attributes, particularly leaf shape and size, crucial to survival and a plant’s adaptation to various climatic challenges [2,3]. As a result, leaves are highly variable among taxa, as well as within each species and even each individual [3,4,5]. While the variability is most often related to functions of photosynthesis [2,6,7] and drought response [8,9,10], leaves have been shown to display variability as a result of numerous climatic influences, including temperature [11,12,13], precipitation [14,15,16] and light conditions [17,18,19,20]. Additionally, the influence of latitude [21,22] and altitude gradient [3] on variations in leaf morphology has been proven. The underlying mechanisms of adaptation are both anatomical and physiological. The structural tissues of leaves, the epidermis, with stomata, mesophyll and vascular tissues, can vary significantly in thickness and shape [23,24]. On the other hand, plant hormones and sensitivity to certain hormone concentrations will influence variable morphological responses [25,26].

High levels of variability, noticeable among different genera as well as within individual species, make leaf morphometric traits valuable in research on population variability [27]. Intraspecific variability, particularly on the population level, might offer beneficial responses to the conditions in the environment, i.e., enable the population to survive [28]. In recent studies, leaf morphometric traits have been used to analyse adaptability of tree species to certain habitats and their specific environmental conditions [29,30,31]. Increased knowledge on the responses of certain population to variability of the environmental conditions can ease management of such populations, possibly aiding conservation of the endangered genera as well [32]. This is particularly important for riparian species, which live in highly diverse environments, where phenotypic plasticity and adaptation are the keys to adapting to ever-changing conditions [33]. One such genus is Salix L., the willows. Willows are shade-intolerant colonisers and true pioneer species, with a rapid growth rate and tendency to vegetative propagation. Furthermore, some species have proven to be morphologically very variable, with pronounced phenotypic responses to environmental factors [34,35,36], such as S. alba L., a common short-rotation forest species, and S. viminalis L., used in remediation efforts [37]. In addition, studies on high-yielding willow genotypes revealed the effect of environment to have greater influence on biomass allocation than the genetic background [38], as well as the influence of specific leaf surface structure on the frost tolerance of the clones [39].

Morphological studies have also been conducted on economically insignificant but ecologically relevant species, such as S. reticulata L. [40] and S. herbacea L. [41], with both species demonstrating a correlation between leaf morphometric traits and the climatic conditions. Similar results were reported for the closely related genus Populus L., the poplars. In an earlier research, black poplar (P. nigra L.) had demonstrated significant phenotypic plasticity of the leaves, providing the species with response to environmental variations [42]. Furthermore, the genetic analysis of P. szechuanica var. tibetica C.K. Schneid. ascertained that certain genes were exclusively expressed in lower altitudes, thus revealing the influence of the environment on the phenotypic expression [43].

A prime example of ecologically significant willow species is the almond-leaved willow (S. triandra L.). Native to Eurasia [44], it grows as a tall shrub, rarely a small tree. The defining feature of the species are almond-shaped leaves with parallel edges, explaining the synonymous scientific name, S. amygdalina L. (from Greek amýgdalon, almond). An important distinctive feature is kidney- or half-heart-shaped stipules, which are permanent and most developed on the top of this year’s shoots and on the suckerings [45]. As with other Salix species, almond-leaved willow is a dioecious with unisexual flowers and is pollinated by insects, with seed dispersal aided by water and wind [46]. Another identifying feature is the male flowers, which are distinguished by having three non-coherent stamens each, giving the species its name (from Latin triandra, three stamens). Growing in a wide variety of habitats, on a wide range of soil types and in various altitudes, almond-leaved willow is empirically adaptable, although previous research to support the observed population variability is unknown to us. To this date, this species has been a part of a research into Salix germplasm [47] and some research into leaf anatomy [48,49] and sprout regeneration [50] has been conducted. Most recently, the pollen grain morphology has been analysed and the species has been separated from the other species within the subgenus Salix, but remained morphologically similar to numerous species from subgenera Chamaetria and Vetrix [51]. Furthermore, although previously divided into numerous subtaxa, the species is now considered unified, with subtaxa genetically classified as synonyms [47] and the analysis of the leaf’s wax layer additionally confirmed these various taxa as a single, unambiguous species [48].

In this study, the phenotypic diversity of almond-leaved willow populations was examined. Plant material, in the form of leaves, was collected from 12 populations growing in diverse climatic conditions. Various statistical methods were employed to analyse leaf morphology variability and test correlations. The hypotheses we tested in this research are as follows: (1) the leaf morphometric traits are positively influenced by favorable environmental conditions; (2) leaf traits are correlated with both geographic location and environmental conditions; (3) intra- and inter-population variability are demonstrated by variability of leaf phenotypic traits; and (4) population structure of almond-leaved willow is under the influence of neutral and adaptive evolutionary processes.

2. Materials and Methods

2.1. Plant Material

Twelve natural populations of almond-leaved willow were selected for the research (Figure 1, Table S1). The habitat conditions varied greatly, with altitudes from 19 to 352 m.a.s.l. and annual precipitation from 723 mm to 1420 mm. Four populations were located in the floodplains of relatively large rivers, growing on alluvial, well-saturated soils. Remaining populations were located in the higher altitudes, with variable distances from the water’s edge and mostly along smaller karstic rivers and rivulets. In each population, five to seven almond-leaved willow shrubs were selected for morphometric analysis. Sampling was conducted during the vegetation period of 2021, after the leaves were fully developed. From each shrub, 20 shoots (Figure 2A), on the external part of the crown and well-lit, were cut and placed in Ziploc plastic bags. Bags were kept in portable freezer containers to additionally prevent them from drying up and deforming. Once collected, samples were taken to the laboratory, and pressed between newspapers and herbarised until the morphometric analysis was conducted.

2.2. Environmental Data

Data on climatic conditions in the area of the studied populations were obtained from the WorldClim2 database with a spatial resolution close to a square kilometre [52]. All 19 bioclimatic variables were included in the study; (Table S1): BIO1 (Annual Mean Temperature); BIO2 (Mean Diurnal Range (Mean of monthly (max temp – min temp)); BIO3 (Isothermality (BIO2/BIO7) (×100)); BIO4 (Temperature Seasonality (standard deviation ×100)); BIO5 (Max Temperature of Warmest Month); BIO6 (Min Temperature of Coldest Month); BIO7 (Temperature Annual Range (BIO5–BIO6)); BIO8 (Mean Temperature of Wettest Quarter); BIO9 (Mean Temperature of Driest Quarter); BIO10 (Mean Temperature of Warmest Quarter); BIO11 (Mean Temperature of Coldest Quarter); BIO12 (Annual Precipitation); BIO13 (Precipitation of Wettest Month); BIO14 (Precipitation of Driest Month); BIO15 (Precipitation Seasonality (Coefficient of Variation)); BIO16 (Precipitation of Wettest Quarter); BIO17 (Precipitation of Driest Quarter); BIO18 (Precipitation of Warmest Quarter); and BIO19 (Precipitation of Coldest Quarter). The bioclimatic variables and altitude were selected to describe the environmental characteristics of the sampling sites for the principal component (PC) analysis and for the calculation of environmental distances. Although some bioclimatic variables were highly correlated, all were included in the further statistical analysis, as they describe environmental characteristics of the studied populations and interpret the differences on the population level to a higher degree.

2.3. Morphometric Analysis

Two leaves from the central part of the shoot were selected at random for further analysis, as those are generally considered to be the most uniform ones [53] with a total of 40 leaves per shrub. In total, 68 individuals were sampled, resulting in morphological data for 2720 leaves. The measurements were carried out using the WinFolia program [54]. Nine phenotypic traits were selected for measurement (Figure 2B). The shape of the leaf was described as the angle closed by the main leaf vein and the line defined by the leaf blade base and a point on the leaf margin, at 10% (LA1) and 25% (LA2) of leaf blade length, and the form coefficient (FC). The remaining six parameters indicated the general leaf dimensions: leaf area (LA); leaf length (LL); maximum leaf width (MLW); leaf length, measured from the leaf base to the point of maximum leaf width (PMLW); leaf blade width at 90% of leaf blade length (LW2) and petiole length (PL). After the above mentioned nine traits were measured, leaves were dried for 24 hours at 105 °C and weighed. The weight of dried leaves was used to calculate specific leaf area (SLA), by dividing mean leaf area by mean leaf dry weight (cm² g⁻¹ dw).

2.4. Statistical Analysis

All data were standardised before analysis to avoid the possible influence of variation resulting from various types of traits. To assess the possibility of conducting multivariate statistical analyses and parametric tests, the symmetry, unimodality and homoscedasticity of data were verified [55]. Assumptions of normality were checked using the Shapiro–Wilk test, and the assumption of homogeneity of variance using Levene’s test.

Arithmetic mean, standard deviation, minimum and maximum value and coefficient of variation were calculated for the particular trait for each population in order to determine the range of their variation. To detect the level of inter- and intra-population variability, hierarchical analysis of variance was used. The analysed factors were populations and shrubs within populations (shrub factor nested inside the population factor). In addition, statistically significant differences between all pairs of populations were identified using Fisher’s LSD multiple comparison test, at p ≤ 0.05. Descriptive statistics and hierarchical analysis of variance were carried out using the STATISTICA software package version 13 [56].

Specific leaf area (SLA) for all of the 12 studied Salix triandra populations was shown in a form of a bar chart, with additional vertical error bars representing standard deviations. One-way analysis of variance was used to test differences between populations.

To analyse structure of the studied populations, several multivariate statistical methods were performed. The K-means method was applied to detect phenotypic structure and define the number of K-groups that best explained the morphological variation of populations [57,58,59,60]. If the proportion of a specific population was equal to or higher than 0.7, that population was assumed to belong to one cluster, and if it was lower than 0.7, that population was considered to be of mixed origin [60]. A dendrogram of the closest Euclidean distances on the basis of the unweighted pair-group method using arithmetic means (UPGMA) was constructed to check the structure between the studied populations. Furthermore, discriminant analysis was performed to evaluate the utility and importance of measured leaf traits by determining which were most useful in maximally discriminating the populations, and to eliminate possible redundant variables. The proportion of individuals correctly classified into the studied populations was determined using classificatory discriminant analysis. The above multivariate statistical analyses were conducted using the “MorphoTools” R scripts in R Version 3.2.2 [61] following the manual of Koutecký [62].

To test correlations between morphometric, geographic and environmental data three different matrices were calculated. Climate data [52,63] and altitude values, retrieved from GPS data recorded during fieldwork, were used to calculate the environmental distance matrix. Environmental differences were calculated as the Euclidian distance between the population means for the first three principal components of the principal component (PC) analysis. Squared Mahalanobis distances between the populations were computed to obtain a matrix of morphometric distances among the studied populations. Geographic distances were calculated from the latitude and longitude of the site of sample collection. Finally, to assess isolation by distance (IBD) and isolation by environment (IBE), response matrix (morphological differences) was compared to the two predictor matrices (climate differences and geographic distance) using simple Mantel tests [64,65,66]. The significance level was assessed after 10,000 permutations, as implemented in NTSYS-pc Version 2.21L [67].

3. Results

3.1. Climate Differences among Sampling Sites

Nineteen bioclimatic variables and altitudes were used to describe the environmental differences between 12 studied populations (Table S1). The first principal component (PC1) explained 67.53% of the total variation and clearly separated the karstic (P02, P03, P04, P05, P06, P07 and P08) populations from higher altitudes from other continental (P09, P10, P11 and P12) and sub-Mediterranean (P01) populations (Figure 1B and Table S2). Named principal component demonstrated strong negative correlation with altitude, BIO3, BIO12, BIO13, BIO14, BIO16, BIO17, BIO18 and BIO19 and positive with BIO1, BIO4, BIO5, BIO8, BIO9 and BIO10. The second principal component explained 19.45% of the total variation and was negatively correlated with BIO6 and BIO11, while positively correlating with BIO7.

3.2. Populations’ Phenotypic Diversity

Our results clearly demonstrate high phenotypic diversity of the almond-leaved willow populations (

Table 1. Descriptive statistics for analysed traits of Salix triandra leaves from 12 studied populations. Populations’ acronyms as in Figure 1. Leaf morphometric traits acronyms as in Figure 2. Maximal and minimal values for arithmetic mean (M) and coefficient of variation (CV) are highlighted with red and green colour, respectively. M—arithmetic mean; SD—standard deviation; Min—minimal value; Max—maximal value; CV—coefficient of variation (%).

Population ID	Descriptive Parameter	LA (cm2)	FC	LL (cm)	MLW (cm)	PMLW (cm)	LW2 (cm)	LA1 (°)	LA2 (°)	PL (cm)
P01	M	5.47	0.37	5.75	1.33	2.78	0.46	28.95	20.17	0.92
	CV (%)	45.15	21.20	24.56	21.87	30.05	28.03	20.21	20.78	32.94
P02	M	4.00	0.39	4.85	1.18	2.02	0.34	34.88	22.79	0.76
	CV (%)	23.71	19.04	14.74	14.98	18.51	32.79	9.81	15.15	21.47
P03	M	3.84	0.37	4.91	1.11	2.09	0.33	32.49	21.05	0.56
	CV (%)	26.94	16.78	16.13	15.23	19.86	31.17	11.12	12.56	23.68
P04	M	2.60	0.39	3.90	0.95	1.78	0.31	32.58	20.99	0.38
	CV (%)	21.48	11.04	13.28	11.85	15.38	17.95	9.46	9.70	23.08
P05	M	3.75	0.40	4.60	1.15	2.06	0.34	34.80	22.67	0.72
	CV (%)	24.15	13.13	11.11	16.02	16.80	31.12	9.00	10.60	15.91
P06	M	4.47	0.40	5.00	1.26	2.08	0.38	35.51	23.29	0.93
	CV (%)	29.41	12.84	15.23	17.71	19.60	27.87	10.00	10.49	22.95
P07	M	3.48	0.36	4.69	1.03	2.02	0.32	32.39	20.73	0.83
	CV (%)	25.35	11.45	13.20	16.57	19.84	20.75	9.27	12.87	23.37
P08	M	3.95	0.41	4.63	1.18	1.86	0.38	37.01	24.01	0.64
	CV (%)	23.13	8.99	11.75	13.46	21.91	25.41	7.78	8.15	18.25
P09	M	5.08	0.34	5.78	1.20	2.33	0.36	32.56	20.10	0.84
	CV (%)	32.14	9.42	15.58	16.05	31.06	31.31	10.90	10.45	26.81
P10	M	3.95	0.33	5.25	1.04	2.16	0.32	31.70	19.04	0.80
	CV (%)	21.62	15.03	11.51	14.94	18.60	25.55	12.31	13.87	23.98
P11	M	7.76	0.35	7.00	1.58	2.48	0.47	34.59	22.22	1.21
	CV (%)	25.04	18.32	13.28	17.04	22.99	26.19	10.18	13.50	20.77
P12	M	6.75	0.33	6.90	1.39	2.58	0.39	33.49	19.89	1.02
	CV (%)	23.23	20.04	12.37	17.22	22.98	28.10	17.06	18.65	21.82
Total	M	4.41	0.37	5.15	1.19	2.15	0.36	33.67	21.62	0.78
	SD	1.75	0.06	1.10	0.24	0.55	0.11	4.31	3.23	0.26
	Min	1.47	0.17	2.81	0.58	0.92	0.23	17.00	12.00	0.13
	Max	16.06	0.66	10.47	2.35	6.40	0.89	48.00	35.00	2.18
	CV (%)	39.74	16.78	21.39	20.37	25.51	30.82	12.80	14.92	33.50

). When observing the values of morphometric traits’ analysis, three populations stand out. The karstic population P04 was characterised by having minimum values for six out of the nine measured traits (LA, LL, MLW, PMLW, LW2 and PL). On the other hand, the lowland population P11 demonstrated the highest values for five traits (LA, LL, MLW, LW2 and PL), whereas the three traits describing the leaf shape (FC, LA1 and LA2) had maximum values in the karstic population P08. The minimum values for the leaf shape traits were measured in populations P10 (LA2 and FC), P12 (LA2) and the sub-Mediterranean population, P01 (LA1). Four karstic populations (P02, P03, P07 and P08) and one lowland (P09) had intermediate values for all of the nine parameters. The most variable population was undoubtedly the single sub-Mediterranean population P01, with seven out of nine traits demonstrating the highest values of the coefficient of variation for both traits describing the shape and the size of the leaves (LA, FC, LL, MLW, LA1, LA2 and PL). The least variable population, with minimum variability for four traits (LA, MLW, PMLW and LW2), was P04.

When observing the SLA values (Figure 3), P04 demonstrated lowest mean values (93.82 cm² g⁻¹), whereas P12 demonstrated the highest (162.73 cm² g⁻¹). Mean data distribution followed a pattern of the distribution for LA, with P11 and P01 having the second and third highest values and P07 as the second lowest value (Figure 3). In addition, for values of SLA one-way analysis of variance (ANOVA) showed significant differences between studied populations.

Significant differences between all of the researched populations were determined for all nine morphometric leaf traits (

Table 2. Results of the hierarchical analysis of variance. Leaf morphometric traits acronyms as in Figure 2.

Trait	Components of the Variance	df	F	Percent of Variability	p-Value
LA	Among populations	11	15.28	53.31	<0.01
	Within populations	56	13.75	15.43	<0.01
	Error			31.26
FC	Among populations	11	2.66	14.41	<0.01
	Within populations	56	26.94	41.85	<0.01
	Error			43.74
LL	Among populations	11	19.77	57.53	<0.01
	Within populations	56	12.12	12.26	<0.01
	Error			30.21
MLW	Among populations	11	6.78	38.55	<0.01
	Within populations	56	21.60	26.82	<0.01
	Error			34.63
PMLW	Among populations	11	6.24	24.48	<0.01
	Within populations	56	10.74	21.54	<0.01
	Error			53.98
LW2	Among populations	11	3.77	15.65	<0.01
	Within populations	56	12.06	23.81	<0.01
	Error			60.54
LA1	Among populations	11	3.39	14.75	<0.01
	Within populations	56	25.34	43.36	<0.01
	Error			41.89
LA2	Among populations	11	2.96	16.40	<0.01
	Within populations	56	33.39	46.41	<0.01
	Error			37.19
PL	Among populations	11	11.63	52.82	<0.01
	Within populations	56	24.13	22.61	<0.01
	Error			24.56

). A larger proportion of total variability of the traits related to leaf length and surface area (LA, LL and PL) was assigned to inter-populational variability, whereas the traits describing the leaf shape (FC, LA1 and LA2) were defined by having the larger proportion of total variability linked to intra-populational differences, although a large proportion was linked to error as well. The three traits related to leaf width (MLW, PMLW and LW1) showed no clear trend in distribution of proportion of total variability among the three components of variance. PMLW and LW2 had the largest proportion assigned to the error component, whereas the remaining proportion was similarly divided between intra- and inter-populational variability. MLW, on the other hand, had proportions of total variability almost equally divided between the three components of variability, 38.55%, 26.82% and 34.63% for inter-, intra-populational and error components, respectively.

K-means clustering inferred the origin of the populations and revealed three clusters, marked in green, blue and yellow in Figure 1A. Easternmost populations of P11 and P12 demonstrated origin completely assigned to cluster A (green), whereas population P01 showed a partial cluster affiliation belonging to cluster A (proportion of membership: 0.6) and cluster B (proportion of membership: 0.4). cluster B (blue) was defined by karstic populations of P02 (proportion of membership: 0.9) and P07 (proportion of membership: 1.0), and the lowland populations of P09 (proportion of membership: 0.8) and P10 (proportion of membership: 0.7). Additionally, three other karstic populations: P03, P05 and P06, showed mixed origin. The third cluster (C, yellow) encompassed two karstic populations, P04 (proportion of membership 1.0) and P08 (proportion of membership: 0.7), with P03, P05 and P06 demonstrating a certain proportion of their origin assigned to Cluster C.

As revealed by UPGMA clustering, researched populations grouped into two clusters with Cluster 1 further divided into two subclusters (Figure S1). This division shares great similarities with the Clusters A, B and C, revealed in the K-means method. The first subcluster of Cluster 1 comprised of two lowland populations, P09 and P10, as well as the karstic populations P03 and P07. The second subcluster was defined by having the remaining karstic populations, with P02 and P05 as being the most similar, whereas the P04 demonstrated the greatest difference within this second subcluster. In the Cluster 2, the easternmost populations P11 and P12 grouped with P01, indicating a lack of clear geographic structuring, but potential great influence of common ecological conditions.

To examine which of the leaf traits were statistically significant between the individual populations, a multiple comparison Fisher’s LSD test was performed (Table S3). The unique character of the P01 population was revealed, with this population having the largest number of significantly different leaf trait values compared to all remaining populations. The easternmost populations, P11 and P12, were highly significantly different compared to the other populations. However, it was found that only four traits were significantly different between populations P01 and P12, while P01 and P11 differed significantly in six traits. Furthermore, populations P11 and P12, also had only four significantly different traits, confirming their similarity. The karstic population, P05, shared common characteristics with populations P02 and P03, as they did not show significant trait differences. In addition, P02 and P05 showed similarity with populations P07 and P08, as they did not have differentiating traits. P04 and P08, although belonging to the same cluster and subcluster, were significantly differ one from another in six out of nine measured traits. When observing overall results, petiole length (PL) proved to be the most significantly differentiating leaf trait.

Six out of nine measured leaf traits have shown to be significant in discriminating the researched almond-leaved willow populations. The greatest discriminating power was noted for traits PMLW and PL, with partial Wilks’ lambda values of 0.49 and 0.50, respectively (Table S4). Other traits’ significance, in descending order, is as follows: LA, LA2, LL and FC. Traits MLW and LW2, as well as LA1, have not shown significant discriminant power among the populations. Overall, for nine traits and 12 populations, the canonical discriminant analysis resulted in nine discriminant functions, with the first three functions demonstrating eigenvalues greater than 1. The first discriminant function, accounting for 63.8% of total variability, was the best discriminator between the easternmost populations of P11 and P12 and the remaining 10 populations (Figure 4 and Figure 5). The next two functions contributed significantly less to the overall variability, accounting for 17.8% and 8.4% of total variability, respectively. The second function discriminated the P04 population best, and the third function discriminated the P01 population best. Generated values of the discriminant functions are shown on the shrub level, which are grouped based on the population affiliation. A very clear separation of populations P11 and P12 along the first axis and populations P04 and P01 along the second and third axis was observed. The classification accuracy for all of the populations was 80.9%. The highest classification accuracy was noted for individuals from the population P11 (100%), whereas the accuracy for the individuals of the population P06 amounted to only 50%.

The analysed populations showed significant correlations between both geographic and phenotypic distances (isolation by distance (IBD), (r = 0.516, p = 0.0077) (Figure 6A)) and between environmental and phenotypic distances (isolation by environment (IBE), (r = 0.493, p = 0.0037) (Figure 6B)).

4. Discussion

Almond-leaved willow research has been scarce, with generally little previously published data on the species [47,48,49,50,68]. Krüssmann [68] reported leaf length (LL) values in the range of 5 to 10 cm. Similarly, Bartha [69] reported leaves to be 3–10 (–15) cm long and (2−) 3–7 times narrower than long. Our mean values of LL (5.15 cm) and MLW (1.19 cm) fall within these data ranges. In addition, mean petiole length (PL) measured for the 12 researched populations was 0.78 cm, which also confirms the previously reported values for PL in range of 0.5–1.5 cm [69].

When observing the variability of leaf traits, a difference between traits related to leaf size and leaf shape can be noticed. The range of coefficient of variation values for leaf size related traits is 20.37% to 39.74%, whereas the leaf shape traits are less variable, with a range of 12.80% to 16.78%, i.e., traits describing leaf size are more variable. The same pattern has been confirmed for other Salix species as well [40,41]. Willows are generally a variable genus, with leaf morphology and anatomy showing high levels of adaptability [39]. Almond-leaved willow is no exception, being previously noted as extremely variable in terms of leaf and shoot morphology [45]. Similarly, a number of poplar (Populus spp.) species, also members of the Salicaceae family, have been shown to have a high degree of leaf size variability [42,45,70,71,72,73]. This trend seems to extend to other pioneer species, including black (Alnus glutinosa (L.) Gaertn.) and grey (A. incana (L.) Moench) alder [60,74,75], and paper (Betula papyrifera Marsh.), silver (B. pendula L.) and yellow birch (B. alleghaniensis Britton) [76,77,78]. As a pioneer species, fast establishment in a new environment is the key to survival and may depend on the morphology of the leaves, primarily leaf area [42]. This would make leaf area the key trait, which could explain the high variability of the trait in this research. In addition, crown architecture, with the goal of preventing self-shading, is an important trait of pioneer species, enabling them the optimal growth conditions [79].

When observing individual populations of Salix triandra, P01 (Mirna) was most variable, with the most variable values of both leaf size and leaf shape traits. This could be attributed to the relatively high humidity of the sub-Mediterranean climate, as well as higher temperatures, two factors that positively influence leaf variability [80]. This is confirmed by the bioclimatic data for the given population, as it is the warmest out of all studied populations (13.3 °C) and has a relatively high precipitation rate (1015 mm) (Table S1). In addition, higher variability of this population could be the result of increased nutrient content in the soil, causing more variation in the production of leaf tissue. The effect of increased nutrient content in soils has a stronger effect on shrubs than on trees, with shrub species, such as almond-leaved willow, positively reaching to higher nutrient levels [81]. Unlike previously mentioned sub-Mediterranean population, karstic population P04 (Vitunj) demonstrated lowest variability on the population level, with the second least variable populations being another karstic population, P08 (Zagorska Mrežnica south). Unlike lowland populations, these populations grow on nutrient-poor soils, out of the reach of the floodwaters, and are often geographically more separate, making gene flow less likely. Furthermore, more demanding conditions may have favoured a smaller number of individuals, leading to reduction of the variation in the populations’ gene pool.

Dioecious plants are known to display modest dimorphism in vegetative and phenological traits [82] and such differences might have influenced leaf morphology of almond-leaved willow as well. The gender-based differences are well-documented, but primarily relate to plant physiology [83], with influence on carbon gain [84], stress tolerance [85], photochemical output [86] and herbivory defence [87]. In some cases, gender-related plant responses are influenced by the environment [85,88,89]. In addition, reported data on the influence of gender on morphology is inconsistent and vary from having influence [90,91,92,93], being of uncertain influence [94,95,96] to no influence at all [97]. These inconsistencies can be particularly observed in the Populus genus, closely related genus to willows. Female plants of Chinese white poplar demonstrated lower values for specific leaf area (SLA) and leaf length (LL) [98], with clear distinction between the genders. On the other hand, P. trichocarpa Torr. et Gray ex Hook. and P. balsamifera L. demonstrated no differences among them [99]. In the Salix species, the results are inconclusive, with S. udensis Trautv. et C.A. Mey. not demonstrating any differences in growth [100], whereas S. glauca L. male plant demonstrated better drought tolerance due to lower stomatal conductance and transpiration rates [101]. Considering the high levels of inter-population variability measured in this research, the influence of gender is likely lesser than that of the environment. With the varying levels of influence of gender on the plant morphology in mind, we cannot exclude the possibility of such influence on our populations. Such investigation would be beneficial for broader understanding of the interactions between the genetic background, the environment and gender-related differences of the species in the future.

The populations included in this research have demonstrated inter-population morphometric variability of traits related to leaf size greater than expected, at the same time indicating low intra-population variability. Such great trait variability among small populations can be attributed to uniform habitat conditions in each population, as well as the founder effect, which most likely influenced the phenotypic variability. The founder effect is a specific form of establishment of a new population by a limited number of individuals, inevitably representing only a small fraction of the original population variability [102]. Almond-leaved willow is known for its propensity to vegetative propagation from broken-off branches which, in addition to habitats along the rivers, could explain the ease of establishing a new colony from a singular or small number of individuals. In addition, the size and large number of seeds produced by willows mean some of the seeds could easily travel downstream or be carried by the wind to a favourable location with stable conditions, where they could grow into new individuals. Once sexually mature, these individuals could continue to expand the population not only by vegetative self-propagation, but by seeds as well, thus forming a new population. A small population size growing in stable and uniform environmental conditions additionally guarantees uniform leaf traits among individuals within each individual population. This theory is supported by similar occurrences in other pioneer species, with leaf morphology correlating to the specific environments [78], demonstrating additionally greater inter-population variability [76]. Interaction between genetic and environmental factors often results in variation in leaf size and shape traits, which play a crucial role in a plant’s ability to capture light and photosynthesise, i.e., in overall plant’s fitness [16,27,103,104]. Accordingly, adjustment of leaf morphology by the plants in the variable environmental conditions is the main mechanism in plant adaptation [105]. High levels of variability in these traits, therefore, can be expected and are supported by our findings.

Grouping of the populations into three clusters, with grouping predominantly explained by common ecological conditions, was revealed. The easternmost lowland populations, P11 and P12, and westernmost, mixed-origin population, P01, grouped together, signifying the common denominator to be the ecological conditions, e.g., habitat. P11 and P12 grow in the riparian zone of river Sava, one of the larger rivers of Pannonian catchment, which forms extensive floodplains. Additionally, these populations show high annual mean temperatures of 11.2 and 11.3 °C, respectively. The westernmost population, P01 (Mirna), is named after a river in the northern Adriatic catchment and grows under the influence of the sub-Mediterranean climate, characterised by the highest mean annual temperature, with 13.3 °C. It is therefore possible to presume a river regime of flooding and the high temperatures conditions, in addition to fluvial soils, form the constant and similar conditions, favouring similar morphological adaptations of otherwise disjunct populations. In addition, nutrient-rich soils and constant conditions of soil humidity can directly influence the repeating pattern of trait values in populations from similar habitats, as it has been previously reported in numerous studies for a wide array of living organisms [106]. On the other hand, the karstic populations included in this research did not form a singular, well-defined cluster. These small populations grow on mostly gravelly, often skeletal soils in less stable conditions, with fluctuating water levels between winter and summer. As a result, although geographically close, these populations have adapted to micro-locations and conditions and thus some of them, such as the populations P04 (Vitunj) and P06 (Zagorska Mrežnica south), have significantly different values for most of the measured traits. This is supported by our field observations for population P04 (Vitunj). Unlike other karstic populations, almond-leaved willows in this micro-location grow on a wet meadow, close to but not along the rivulet Vitunjčica. The fact that the plants do not grow directly on the water’s edge most likely led to slight but noticeable differences in soil humidity and perhaps soil type, which in turn have caused specific trait values and variability.

Same differentiation between the karstic and lowland populations can be corrobo-rated by SLA data, which followed the same trend as LA values. High SLA values in low-land populations, characterised by humid and nutrient-rich soils, and lower values in karstic populations found on drier, nutrient-poor soils, are in agreement with previous data on positive correlation between increased humidity and nutrient content and SLA [107,108,109]. SLA is known to be a good predictor of plant growth, as it correlates with the growth parameters of the whole plant [110]. In addition, SLA is a good indicator of the ability of plants to utilise resources available from the environment, with differences between high-SLA and low-SLA species [111]. High-SLA species are usually found in nutrient-rich habitats, where they can afford the shortened life span of leaves. Low-SLA species, however, usually grow on nutrient-poor habitats and thus have more dry matter per leaf, keeping leaves for a number of seasons [111,112]. Therefore, almond-leaved willow would be classified as high-SLA species, a fact supported by populations found on the nutrient-rich soils having the highest SLA values.

The results of this research indicated the existence of both isolation by distance (IBD) and isolation by environment (IBE) patterns. These models refer to morphological variability resulting from geographical (IBD) or environmental (IBE) distances between the populations [113,114]. In other words, the differences in altitude and climatic conditions, in addition to geographical distances between the populations, were great enough to bring about variations in phenotypic expression of almond-leaved willow shrubs, with a trend of similarities between ecologically and geographically closer populations. The certain gene flow pattern can influence species’ populations in various ways, from increasing the genetic variation and population size, to weakening adaptive genetic combinations, leading to a reduction of population size [115]. In small populations, such as those covered in our study, gene flow between similar environments may increase population size and introduce new, locally beneficial alleles [116,117], which will allow plants to adapt to rapidly changing habitats due to climate change. The influence of both patterns on phenotypic variability of woody species has been previously reported by several authors, where geographically closer and environmentally similar populations had similar morphotypes [60,118,119,120,121,122]. This was further confirmed in a comprehensive review by Sexton et al. [115], who discovered that both models have a significant influence in 60% of the cases. This means that gene flow between populations, and thus their phenotypic similarity, is greater at shorter distances and analogous environments. Therefore, IBD in our study might be caused by limitations in seed and pollen dispersal and drift, while IBE could be present due to non-random mating originating from differences in timing of phenological processes and local adaptation caused by strong selection [115].

5. Conclusions

The traits related to leaf size were more significant in describing differences between populations and also had the highest variability. The large differences among the populations are the result of both geographical and environmental distances between them. In addition, considering the small size of the populations, the high variability among the populations can be attributed to the founder effect, which is very likely considering the ease of the almond-leaved willow’s propagation. The results of the multivariate statistical methods revealed three groups of populations defined by common environmental conditions and geographical proximity.