Genetic, geographic, and climatic factors jointly shape leaf morphology of an alpine oak, Quercus aquifolioides Rehder & E.H. Wilson

Abstract

Key message

Leaf symmetry and leaf size are explained by genetic variation between and within lineages and to a lesser extent by climatic factors, while leaf asymmetry can only be partly explained by geographic factors in Quercus aquifolioides Rehder & E.H. Wilson.

Context

Leaves are the primary photosynthetic organs of plants, and their morphology affects various crucial physiological processes potentially linked to fitness.

Aims

We explored the variation in leaf morphology of an alpine oak, Quercus aquifolioides, in order to examine its relationship to genetic, geographic, and climatic factors.

Methods

We conducted a genetic survey using 25 nuclear microsatellites. Based on Bayesian clustering analysis, 273 sampled trees from 29 populations of Q. aquifolioides were assigned to two lineages that correspond to the Western Sichuan Plateau-Hengduan Mountains (WSP-HDM) and Tibet geographic areas, with some individuals showing mixed ancestry. To undertake morphological analyses, we collected 1435 leaves from these trees and characterized them in terms of 13 landmarks. The metric dimensions of these leaves were digitally captured in the two-dimensional coordinates of these landmarks, then divided into leaf size and symmetric and asymmetric components of leaf shape. To analyze how different components of leaf morphology vary across lineages, we employed Procrustes Analysis of Variance (ANOVA), two-block partial least-square analysis (2B-PLS), and several other multivariate analysis approaches. We also applied distance-based redundancy analysis (dbRDAs) to explore relations between leaf morphology and genetic, geographic, and climatic factors.

Results

Multivariate analysis indicated significant differentiation in leaf symmetric shape components and leaf size between the WSP-HDM and Tibet lineages, while the mixed individuals were morphologically intermediate. The dbRDA analysis showed that most of the variation in symmetric components and leaf size was explained by genotypic effects, with the symmetric components of leaf shape being also significantly explained by geography and climate; however, variation in asymmetric components is only very weakly explained by geography.

Conclusion

Our results demonstrated that leaf morphological variation in shape and size across Q. aquifolioides geographic range is related to both its genetic differentiation and to a lesser extent to climatic factors. We discuss how these patterns could be interpreted in terms of both geographical isolations among and within lineages, and possible adaptive responses for particular traits, in contrast to asymmetric variation.

1 Introduction

As the photosynthetic and transpirational organs, leaves play a critical role in the physiological process of plants through exchanging air and maximizing light capture (Nicotra et al. 2011; Ferris 2019). In addition, leaf trait variation (including leaf size, shape, integrated traits, and allometric relationships) might reflect plant adaptations and plastic responses to different environments (Viscosi 2015; Cuevas-Reyes et al. 2018; Ferris 2019; Klápště et al. 2020; Maya-García et al. 2020). Several studies, largely focused on model species, have shown that leaf variation is strongly influenced by both genetic and environmental factors (Barkoulas et al. 2007; Ferris et al. 2015; Fritz et al. 2018). The most commonly reported is that leaf morphology is determined by multiple genes simultaneously (Klingenberg 2010; Tian et al. 2011; Chitwood et al. 2013; Baker et al. 2015; Fu et al. 2017). Some studies show that geographical isolation may initiate leaf morphological differentiation (Krauze-Michalska and Boratyńska 2013; Maya-García et al. 2020). Other studies suggest leaf traits vary with climatic factors (Royer et al. 2008; Fu et al. 2017; Guerin et al. 2012). However, testing jointly different factors for a better understanding of their relative effects on leaf morphological variation is rarely performed. This could provide new insights into processes of divergence and adaptation within species, and thus help developing conservation strategies under climate change.

Both traditional morphological measurements and modern geometric morphometric methods (GMMs) have been widely applied in recent leaf morphological analyses. In traditional morphometric methods, measured lengths and widths of leaf structures and distances between certain landmarks (i.e., points that can be located precisely) are typically subjected to multivariate statistical analyses (Rohlf and Marcus 1993), but the original geometric relationships among the selected landmarks may not be reconstructable (Mitteroecker and Gunz 2009). GMMs are based on Cartesian landmark coordinates and analyses the relative positions of morphological landmarks to represent each specimen (Klingenberg 2011). Thus, results of GMM-based statistical analyses preserve distances between shapes, and types of graphs (e.g., transformation grids) of GMMs can be used to visualize the trend of shape changes (Mitteroecker and Gunz 2009). Furthermore, GMMs can be applied to distinguish symmetry (the repetition of parts in different positions and orientations to each other), bilateral asymmetry (the differences between corresponding parts on the left and right sides), and allometric relationships, which reflect the covariation of size and shape, and are also critical components of morphological studies (Klingenberg et al. 2002; Klingenberg 2011, 2015; Viscosi et al. 2012). Separate analyses of these components of leaf variation can help to quantifying the symmetric components of shape traits, which could be more heritable than the asymmetric components of bilateral structures likely caused by developmental instability (see Viscosi 2015 for GMMs applications to white oaks). Thus, GMMs have been increasingly applied in the analysis of morphological variation (e.gAlbarrán-Lara et al. 2010; Viscosi 2015; Liu et al. 2018; Tucić et al. 2018; Jiang et al. 2019).

Oaks (Quercus L., Fagaceae) are mainly located in temperate forests in Europe, North America, and Asia (Manos et al. 1999) and 35 species are widely distributed in China (Huang et al. 1999). The morphological discrimination among Quercus species is challenging, because of individuals with intermediate traits in many species, as well as both intra and inter-specific morphological variations partly caused by extensive hybridisation and introgression (Whittemore and Schaal 1991; Bruschi et al. 2003a; Lepais et al. 2009; Leroy et al. 2017; Lang et al. 2018; but see also Abadie et al. 2012 and Lepais et al. 2013 documenting strong reproductive barriers in European white oaks). Leaf morphological traits have been traditionally applied to discriminate oak species in mixed stands in both Europe and North America, either using morphological features alone (e.g., Kremer et al. 2002; Ponton et al. 2004; González-Rodríguez and Oyama 2005; Kelleher et al. 2005), or combining genetic and morphometric approaches to study the morphological differences among species (Gügerli et al. 2007; Viscosi et al. 2012; Beatty et al. 2016; Porth et al. 2016; Rellstab et al. 2016), or how leaf shapes respond to the environmental variation (Maya-García et al. 2020; Ramirez et al. 2020). In contrast, published analyses of leaf shapes of Quercus species with Asian distributions are scarce: Sun et al. (2016) studied leaf morphological variation and its response to high altitude in an evergreen oak in Qinghai-Tibet Plateau. Only one GMMs study in oaks in China by Liu et al. (2018) demonstrated that two congeneric deciduous oak species in sympatric areas were efficiently discriminated, and their variations in leaf morphology correlated with environmental factors.

Quercus aquifolioides Rehder & E.H. Wilson is an alpine species endemic to the Himalaya–Hengduan Mountains biodiversity hotspot. The species is widely distributed across the Hengduan Mountains, southwestern China (Zhou 1993), and plays essential roles in preventing soil erosion, regulating climate, and maintaining ecological stability in this region (Xu and Guan 1992; Zhou and Guan 1992). It occupies diverse habitats at altitudes ranging from 2000 to 4500 m, showing immense adaptability to cold and dry environments (Li et al. 2006). Phylogeography study based on nuclear microsatellite (nSSRs) suggested relatively high genetic differentiation (F_ST ~ 7%) between two regions across the Q. aquifolioides geographic range, which we referred to as the Western Sichuan Plateau-Hengduan Mountains (WSP-HDM) lineage and the Tibet lineage (Du et al. 2017). Furthermore, landscape genomics analysis based on 65 drought-stress related candidate genes revealed contrasting patterns of genetic differentiation among populations within both lineages, with stronger evidence for local adaptation to changing environments in the WSP-HDM, while both Tibet and WSP-HDM lineages showed significant patterns of isolation by distance (Du et al. 2020). The purpose of this study is first to examine whether morphological variation observed within Q. aquifolioides can be significantly partitioned among the WSP-HDM and Tibet lineages, and if so, to assess the relative importance of geographic distribution, neutral genetic processes, and local climatic factors in shaping leaf morphological variation or explaining various components of leaf morphological variation. First, we confirmed the genetic assignment of all sampled individuals to the main geographic lineages using additional nuclear SSR markers and multilocus genotypes. We then used GMMs to analyze different components of leaf morphology (including leaf size, symmetric and asymmetric components of leaf shape) between and within lineages. Finally, we analyzed the relationships between leaf morphological traits and genetic, geographic, and climatic factors using distance-based redundancy analyses (dbRDAs), and discussed how the observed patterns could be linked to the species past evolutionary history.

2 Materials and methods

2.1 Sampling strategy

We sampled 273 individuals in 29 Q. aquifolioides populations throughout the species’ distribution (Fig. 1a). Each population was at least 30 km apart, and all sampled individuals were located at least 5 m from each other to minimize risks of selecting clones. We collected six intact, mature leaves from each individual, one of which was dried and stored in silica gel for DNA extraction, while the other five were used for morphological analyses. We measured the latitude, longitude, and altitude of each population using an Extrex GIS monitor (Garmin) (Appendix Table 4). All sampled specimens were deposited in archives of the Molecular Ecology Laboratory of Beijing Forestry University, China.

Fig. 1

Geographical distribution, locations of study sites and mapping of the Bayesian genetic clusters in Quercus aquifolioides. a Species range (dashed line) and sampling locations (black dots). b The magnitude of Delta K and LnP (K) as functions of K suggesting the existence of two major clusters as the most likely scenario. c Histogram of individual assignments for K = 2 and K = 3. The major lineages are indicated on top and population codes below the histogram. d Geographical distributions of two primary lineages and their mixed individuals based on the clustering analysis (K = 2). Threshold Q values 0.8/0.2 were chosen to assign individuals to lineages or their mixed individuals

2.2 Genetic assignment of the specimens

In total, 25 polymorphic, neutral nuclear microsatellite loci were used for genotyping all sampled trees after initial tests with 85 loci. Among them, 11 loci (MSQ13, QpZAG9, QpZAG16, QpZAG110, QrZAG7, QrZAG11, QrZAG30, QrZAG87, QrZAG96, QrZAG112, and PIE271) were extracted from a previous study (Du et al. 2017) and the other 14 (QmC02241, CN725667, CR627959, GOT011, GOT012, GOT021, GOT040, PIE163, FIR026, WAG066, WAG068, POR017, POR025, and FIR015) were genotyped in this study (Ueno et al. 2008; Ueno and Tsumura 2008; Durand et al. 2010). Detailed information on the primers, amplification, and genotyping procedures is provided in Appendix Table 5.

Bayesian cluster analysis was performed for SSR-genotyped individual with STRUCTURE V2.3 (Pritchard et al. 2000). The program was run for one to 10 clusters (K), with 20 repetitions for each K of 200,000 Markov Chain Monte Carlo cycles (MCMC) following 100,000 burn-in cycles. The most likely number of clusters (K) was determined using ΔK and LnP(K) statistics (Fig. 1b), according to Evanno et al. (2005) and Janes et al. (2017). In order to further explore possible genetic clusters, we provide STRUCTURE plots for different K for visual comparison in Fig. 1c. Threshold Q values 0.8/0.2 were chosen to best combine efficiency and accuracy when assigning individuals to lineages (as defined by Vähä and Primmer 2006). We interpreted the individuals with Q values between 0.2 and 0.8 as representing “mixed” lineage (Fig. 1d).

2.3 Landmark configuration

Leaves were pressed, dried, and scanned with the abaxial surface uppermost using a CanoScan 5600 F scanner (Canon Inc., Japan) at a resolution of 600 dpi. Scanned images were characterized by 13 landmarks (LM) (Fig. 2) on the left and right sides of each leaf, as suggested for other oak species (Viscosi 2015; Liu et al. 2018). Three of these landmarks, LM1-LM3, are unpaired and distributed along the middle axis of leaves, while the others (LM4-LM13) are in pairs and occur symmetrically on both sides of the leaves (Savriama and Klingenberg 2011).

Fig. 2

Configuration of Quercus aquifolioides leaves, showing locations of the 13 features used as landmarks (LM) in this study, with descriptions of the landmarks on the right

Cartesian x and y coordinates of the 13 LMs for each leaf (1,435 leaves in total) were recorded and stored in.txt file format using the Image J program (Abràmoff et al. 2004). To estimate digitizing error, a subsample of 290 leaves (from two randomly selected trees from the 10 sampled in each population) was repeatedly marked a month later, as suggested by Viscosi (2015).

2.4 Morphological analysis of leaf variation

To analyze leaf morphology, the Morpho J software (Klingenberg 2011) was used. Briefly, the landmark data were imported into the software; then, shape and size information was separated by Procrustes superimposition, a method fitting multiple configurations onto a common consensus and removing the components that are not part of shape from the coordinate data (Dryden and Mardia 1998; Klingenberg 2010). In this process, the leaves of each tree were processed to create a mirror image by inverting the sign of the x coordinate of each landmark; then, an average configuration was obtained from the original and mirrored configurations to separate the leaf shape variation into symmetric and asymmetric components (Appendix Fig. 5) (Mardia et al. 2000). Here, we defined symmetric variation as the variation in the averages of original and mirrored configurations, and asymmetric variation as the differences between the original and mirrored configurations (Klingenberg et al. 2002). Next, as commonly recommended, we removed four outliers that deviated sharply from the averages, and created a wireframe resembling a stylised leaf by drawing lines between pairs of landmarks to assist visualisations. Finally, the covariance matrix of each hierarchical level (leaves, trees, and lineages) was created for subsequent multivariate statistical analyses by the Morpho J software. All the original leaf morphological datasets were deposited to figshare repository (Du 2021).

We used Procrustes ANOVA to test relative amounts of asymmetric variation and measurement error among leaves, trees, and lineages (Klingenberg 2003). Variance in size or shape was partitioned into main (lineages), random (trees, leaves), and digitizing error (replicates) components. In addition, the relative amounts of the shape variation were also decomposed into components of the directional asymmetry (reflections) and the fluctuating asymmetry (lineages-by-reflections, trees-by-reflections, and leaves-by-reflections interaction), representing systematic differences and small random deviations between the left and right sides of the blade, respectively (Savriama and Klingenberg 2011). The magnitude of the effects was measured in terms of F ratios and percentages of variance explained by each effect.

Symmetric and asymmetric components of all leaves were subjected to Principal Component Analyses (PCAs) to identify differences in leaf shape among lineages. The significance of differences between lineages was tested by Scheffe post hoc comparisons using the R package agricolae (Mendiburu 2020). Allometric patterns of covariation between leaf size and shape, which are critical components of morphological studies, were assessed by two-block partial least-squares (2B-PLS) analysis (Rohlf and Corti 2000). The covariation was measured by the RV (squared Pearson correlation) coefficient, the significance of which was assessed with 10,000 permutations. In addition, canonical variate analysis (CVA) and discriminant analysis (DA) were used to detect tree-level differences between groups (here, lineages). CVA maximizes the separation of specified groups by ordination analysis with permutation tests using Mahalanobis and Procrustes distances, while DA provides reliable information on differences among groups by a cross-validation procedure with the T² statistic (here, P values for tests with 1000 permutations < 0.0001; Klingenberg 2011).

2.5 Effects of geographic, genetic, and climatic factors on leaf morphology

Next, we investigated the amounts of variation among all populations in leaf morphology (size, and both symmetric and asymmetric components of leaf shape) explained by geographic, genetic, and climatic factors. For this, we applied a series of marginal (full) and conditional (partial) dbRDAs with variance partitioning based on a matrix of Euclidean distances (Legendre and Anderson 1999). All dbRDAs were performed using the capscale function in the R package vegan (Oksanen et al. 2013). The marginal dbRDAs models included all variables, while the conditional models accounted for variation in the selected variables. The first three principal components (PCs) of SSR alleles based on PCA were regarded as representative of neutral genetic variables, which explained more than 60% of the total variance. To acquire relevant climatic variables, we downloaded climatic data from the WORLDCLIM database (http://www.worldclim.org) and mapped the GPS coordinates of sampling sites to 30 arc-second ESRI® climate data grids using ArcMAP (ESRI 2014). The highly correlated climatic variables with a high variance inflation factor (VIF > 0.7) were removed from 19 bioclimatic variables using the R package usdm (Naimi 2013). Four climatic variables (Isothermality (Bio03), precipitation of the driest quarter (Bio17), December precipitation (Prec12), and March wind speed (Wind03)) (Appendix Table 6) and three geographic factors (longitude, latitude, and altitude) were retained as explanatory variables. We first used variance partitioning to construct a marginal model including geographic, genetic, and climatic explanatory variables identified in the forward procedure (geographic, genetics, and climatic factors). As one of our main objectives was to test the relative importance of geographic, genetic processes, and local climate in predicting leaf morphological variation, we then ran three conditional dbRDAs models to test the pure geographic, genetic, and climatic effects. These were conditional models of geographic variables controlling genetic and climatic effects (designated geography + conditional (genetics + climate)), genetic variables controlling geographic and climatic effects (designated genetics + conditional (geography + climate)), and climatic variables controlling geographic and genetic effects (designated Climate + conditional (geography + genetics)). The significance of each independent factor included in the dbRDAs was assessed by the anova function in the vegan package, which computes pseudo-F ratios, variance components, and P values (Oksanen et al. 2013).

3 Results

3.1 A prioriassignment of the individuals to the two lineages

The results of the STRUCTURE analysis showed that the estimated log probability of data increased sharply from K = 1 (LnP(K) =  − 25,600) to K = 2 (LnP(K) =  − 24,700) and then increased slightly from K = 2 to K = 3 (LnP(K) =  − 24,500; Fig. 1b). According to Pritchard et al. (2000) and ΔK (Fig. 1b), this pattern should be interpreted as K = 2 being the most likely number of clusters. This was consistent with STRUCTURE clustering results of Du et al. (2017). The genetic assignments of individuals using K = 3 are similar to those with K = 2 for most populations, with a third possible sub-cluster including three populations located at a short geographic distance from one another (Fig. 1c). However, the main genetic clusters as identified from this Bayesian analysis remain those linked to the Western Sichuan Plateau-Hengduan Mountains (WSP-HDM) and Tibet lineages (Fig. 1c). Twenty-four individuals, mainly members of the BMZ and LZD populations, were genetically admixed from these groups, and thus called “mixed” (Fig. 1d).

3.2 Procrustes ANOVA of leaf size and shape

We subjected the leaf size and shape data to Procrustes ANOVA, hierarchically partitioning variance into leaves, trees, and lineages (Table 1). The digitizing error explained a negligible part of the total variance for both leaf size and shape. For leaf size, most of the total variance was associated with between-trees and among-leaves effects, which explained 47.4% and 40.0% of the total variance, respectively, while lineages explained 12.5% of the total variance. For leaf shape, leaf and tree effects also explained a large part (c. 73.2%) of the total variance, while the differences in leaf shape among lineages were significant but only explained c. 2.7% of the variance. In addition, directional asymmetry (the reflections effect) was insignificant while fluctuating asymmetry, the (lineages + trees + leaves) × reflections effect, was highly significant, explaining c. 23.8% of leaf shape variance.

Table 1 Results of the Procrustes ANOVA of leaf size and shape using all samples

3.3 Leaf morphological variation and allometric patterns

In PCA score plots of the symmetric component, the WSP-HDM and Tibet lineages were partially separated (Appendix Fig. 6), and mixed individuals were scattered between the two lineages, but the lineages (and mixed specimens) overlapped almost completely in score plots of the asymmetric components (Appendix Fig. 6). Similarly, post hoc analysis revealed remarkable differences in the symmetric components and leaf size between the WSP-HDM and Tibet lineages, but no distinctions between groups in the asymmetric components (Table 2). In addition, the associations between variations in shape and size indicated that there were significant allometric patterns in the symmetric but not asymmetric variation (Fig. 3). Graphical reconstruction of the symmetric components (Fig. 3a) showed that as leaf size increased, the leaf shape changed from suborbicular to subelliptical, and the relative length of the petiole declined.

Table 2 Differences in symmetric component, asymmetric component, and leaf size among lineages determined by the Scheffe post hoc comparison test (P < 0.05)

Fig. 3

Scatter plots obtained from two-block partial least squares (2B-PLS) analysis, with 95% confidence ellipses, illustrating the relationship between size (log centroid size) and symmetric component (a) or asymmetric component (b) of the leaf shape of Quercus aquifolioides WSP-HDM, Tibet lineages, and mixed individuals. The thin-plate spline deformation grids in a represent leaf shapes reflecting the negative (− shape) and positive (+ shape) extremes of the PLS axis for symmetric component

3.4 Leaf morphological discrimination

Discrimination among groups also showed significant morphological differences between Tibet and WSP-HDM lineages (Fig. 4; Appendix Fig. 7a, T² = 126.7684; p < 0.0001), with percentages of correctly classified trees in each of their lineages being 80% and 87%, respectively. CVA revealed that Tibet and WSP-HDM lineages formed two groups despite a relatively large overlap among groups, and mixed individuals were scattered between the two lineages along CV axis 1 (CV1, 89% of the total variance) (Fig. 4). However, the two lineages largely overlapped along CV2 (11% of the total variance). Transformation grids showed that leaves of trees of the Tibet and WSP-HDM lineages had subelliptical and suborbicular contours, respectively (and the former were blunter in both apical and basal regions), while leaves of the mixed individuals had intermediate positions along CV1 (Fig. 4). The pairwise comparisons revealed greater differences in leaf shape between the two lineages than between either of the lineages and mixed individuals.

Fig. 4

Scatter plot obtained by canonical variate analysis (CVA) computed on the symmetric component with 95% confidence ellipses. Transformation grids represent average leaf shapes of each lineage

3.5 Effects of geographic, genetic, and climatic factors on leaf morphology

Marginal (full) dbRDA (Table 3 and Appendix Table 7) showed that geographic (longitude and altitude), genetic (genotypic PC2), and climatic factors (Wind03) significantly explained the symmetric components of leaf shape (c. 3.8%, 19%, and 7.1% of total variance, respectively). A small proportion of the variation in the asymmetric components of leaf shape can be significantly explained (c. 2.2%) by geography (latitude and longitude) only, and leaf size was significantly explained by both genetic (genotypic PC2) and local climatic (Wind03 and Bio17) factors (c. 22.5% and 23.3% of total variance, respectively). These results were confirmed by partial RDA, in which the other two predictors were conditioned when testing the significance of the proportions of variance in leaf shape and size explained by every single predictor. For the symmetric component, the variation was significantly explained by genetic (PC2; c. 19%), and to a lesser extent geographic (longitude and altitude; c. 2%) and climatic factors (Wind03; c. 3.5%), while for the asymmetric component, only geography (latitude and longitude) was statistically significant, but explaining little of the total variance (1.8%). The main difference between marginal and partial dbRDAs was that climatic variables (mainly Wind03) significantly contributed to the explanation of leaf size according to marginal but not partial dbRDA (only genetic (PC2) 22.5%; Table 3 and Appendix Table 7). Results of the Mantel test (Appendix Table 8) showed a significant linear relationship between climatic and genetic distances, indicating that climatic effects on leaf size might be largely due to interactions between genetic and climatic factors, rather than climate per se.

Table 3 Results of tests of the association of symmetric component, asymmetric component and leaf size (all samples) with three sets of predictor variables using distance-based redundancy analysis (dbRDA): marginal tests on the left and partial (conditional) tests in which we tested the significance of each variable set, controlling the other two variable sets, on the right

4 Discussion

It is widely accepted that various geographic, genetic, and climatic factors can contribute to changes in leaf morphology (Ferris et al. 2015; Sun et al. 2016; Fritz et al. 2018; Maya-García et al. 2020). In this study, we first demonstrated that Q. aquifolioides lineages identified from their multilocus genotypes differed morphologically, and that mixed individuals among genetic lineages were morphologically intermediate. Using dbRDAs, the relative importance of three types of factors (i.e., geographic, neutral genetic variation, and climatic) was tested on the different components of leaf morphological variability. We showed that all factors were correlated to the variation of different components of the species’ leaf morphology throughout its range. Most interestingly, the variation of symmetric components and leaf size were mainly explained by genotypic variation and to a lesser extent by climatic factors, while the asymmetric components were only weakly related to geographic factors. We discuss these results in the context of the ecological and genetic literature on variation in leaf morphology.

4.1 Possible impact of lineages, geographical isolation on leaf size, and shape variation

Variation in leaf morphological traits is often important among plants of the same taxa (including trees) and even within plant species (Baranski 1975; Blue and Jensen 1988; Bruschi et al. 2003b). In Q. aquifolioides, large and similar proportions of variation (36 to 47%) were observed for leaf shape and size variability among trees and between tree effects explaining slightly more variation than within tree effects for leaf size (Table 1). The same result in white oaks was interpreted as being due both to the fact that leaves from the same trees share the same genes, but also a more homogeneous micro-environment compared to leaves from different trees (Viscosi and Cardini 2011). Between the WSP-HDM and the Tibet lineages, we also detected significant differences in leaf size, as well as relatively small but significant differences in leaf shape (12.5% and 2.7% of variance explained, respectively, in Table 1). This could result from geographical isolation illustrated by a relatively high neutral genetic differentiation (Du et al. 2017, overall F_ST of ~ 7%) that could be related to the geological past history of these mountainous regions. Indeed, Favre et al. (2015) demonstrated that parallel north–south oriented valleys surrounded by high peaks were formed in south-eastern Tibet and north-western Yunnan as a consequence of the mountain uplifts, promoting geographical isolation among populations that could have driven morphological differentiation (Krauze-Michalska and Boratyńska 2013). This interpretation is consistent with the significant proportion of variation explained by the genotypic factor for both the symmetric components of leaf shape and leaf size traits (~ 20% in Table 3, using a more complete set of nSSRs), when using dbRDA to jointly test geographic, genetic, and climatic factors on morphological variation. Besides, Du et al. (2020) also showed strong patterns of isolation by distance overall and within lineages, using a few hundred SNP markers, consistently with genetic drift and limited gene flow between distant populations.

4.2 Impact of environmental pressures on leaf shape and size traits

In the WSP-HDM group of populations, divergent selection pressures were also recently suggested to interpret molecular signals of differentiation at particular drought-stress related candidate genes (Du et al. 2020). This could be linked to the observation that climatic variables explained a significant part, although relatively low, of the symmetric components of leaf shape variation across the species range (Table 3, 7.1% and 3.5% of total variance for marginal and conditional dbRDA tests respectively).

Allometric patterns could indeed reflect plant adaptation to heterogeneous environments, accounting for a large portion of the total variation in plants’ geometric morphology (Klingenberg 1997a, b; Klingenberg 1998; Viscosi et al. 2012), for example, along altitudinal gradients (Bresson et al 2011), where variation in temperature and precipitation occur. Here the climatic explanatory variables included temperature daily oscillations (Isothermality), interaction of precipitation levels with temperatures (during driest quarter and coldest month), and wind speed. Previous research showed that geographical factors may influence species’ allometric relationships (López-Serrano,2005). In this study, we found that although both lineages of Q. aquifolioides had very similar allometric patterns, the extent of shape changes was higher across the WSP-HDM lineage (Fig. 3). This might reflect possible higher environmental heterogeneity and local adaptation evidence in this particular lineage, while comparable allometric patterns were observed in other species of deciduous oaks (Viscosi 2015; Liu et al. 2018). In contrast, we found evidence for an isometric relationship (i.e., non-significant links between size and shape) between the asymmetric components of leaf shape and size in Q. aquifolioides. More studies in Asian Fagaceae with widespread distributions would help determining the extent to which allometric and isometric elements of the leaf shape/size relationship are conserved across different habitats.

Despite the strong similarity of allometric patterns, clear morphological differentiation was detected between the WSP-HDM and Tibet lineages as mentioned previously, mainly related to their leaves’ basal and apical regions. Mixed individuals exhibited an intermediate leaf shape, but pairwise comparisons showed that they also overlapped with the distributions of individuals from the two principal genetic clusters (the WSP-HDM and Tibet lineages).

4.3 Origin of fluctuating asymmetry?

In contrast to allometric patterns of morphological variation, fluctuating asymmetry is generally considered as representative of random phenotypic responses to environmental perturbations, and it has been widely used in evolutionary biology studies as a measure of developmental instability (Møller and Swaddle 1997; Alibert and Auffray 2003; Polak 2003; Graham et al. 2010; Klingenberg 2010). Although fluctuating asymmetry represents a reasonable part of total leaf shape variation in Q. aquifolioides (Table 1, ~ 24% of variance explained), redundancy analyses show that its variation was significantly but weakly related to geographic factors only (marginal dbRDA test: c. 2.2% of total variance; conditional dbRDA tests: c. 1.8% of total variance; Table 3). These results corroborate previous indications that genetic and climatic factors weakly affect asymmetrical variation (e.g., Furuta et al. 1995; Iwata et al. 2002). The absence of more systematic genetic or climatic factors explaining fluctuating asymmetry might indicate that this type of variation occurs in a similar way across Q. aquifolioides populations. It could be due to responses of tree leaves to generally stressful environmental conditions encountered across all populations studied, with habitats characterized by steep high solar radiation slopes of rugged mountains and high elevation (Huang et al. 1999; Tang 2006).

5 Conclusions

Our study of Q. aquifolioides combined morphological, geographic, genetic (nSSRs), and climatic analyses in efforts to obtain an integrated understanding of variation in leaf morphology. Using GMMs analyses to characterize leaf size and shape variation on one hand, and building from previous knowledge on genetic structure among geographical lineages, we observed significant intraspecific variation in Q. aquifolioides leaf shape and size traits across lineages and individual trees. Specifically, we found that leaf symmetry and leaf size could be related to genotypic variation, consistently with lineages differentiation, and patterns of isolation by distance within lineages. In addition, part of the variation in leaf symmetric components could also be interpreted as possible adaptive responses to heterogeneous habitats. In contrast, leaf asymmetry is largely unexplained by genetic and climatic factors, only weakly regulated by geographic factors, consistently to its generally acknowledged developmental instability origin. This calls for caution in using leaf traits for taxonomic purposes: symmetric components and leaf size were suggested for taxonomy but not asymmetric components. Our results provide potentially valuable information for understanding the ecology and evolution of Q. aquifolioides populations. Further studies including controlled experiments could however be useful for testing various hypotheses that could link variable leaf shape characteristics to putative adaptive functions to different environments.

Appendix

Table 4 Geographical location of each population and numbers of individuals within it assigned to indicated lineages based on genotype assignment (Q = 0.8)

Table 5 Information regarding the 25 nuclear simple sequence repeats (nSSR) loci primers used for genotyping

Table 6 Climate variables for sampling localities for Quercus aquifolioides. Bio17 precipitation of driest quarter, Bio03 isothermality, Prec12 December precipitation, Wind03 March wind speed

Table 7 Results of tests of the association of symmetric component, asymmetric component and leaf size on the total sample with several individual predictor variables, using the dbRDA multivariate F-statistic. %VAR, percentage of variance explained by each variable; F, F values; ***P < 0.001; **P < 0.01; *P < 0.05

Table 8 Results of Mantel and partial Mantel tests of pairwise relations between genetic distance (F_ST/(1 − F_ST)), symmetric component, asymmetric component, and leaf size with geographic and environmental distances, using all samples. ***P < 0.001; **P < 0.01; *P < 0.05

Fig. 5

Results of generalized Procrustes analysis of the leaf shape of Quercus aquifolioides based on Cartesian x and y coordinates of 13 landmarks (LMs): a using the full raw coordinate matrix; b and c using the separated symmetric and asymmetric components, respectively

Fig. 6

Results of Principal Component Analysis (PCA) of the symmetric component (a) and asymmetric component (b) of leaves of two Quercus aquifolioides lineages (WSP-HDM and Tibet) and mixed individuals. Scatter plots of PC1 and PC2 scores, with 95% confidence ellipses in a. Transformation grids for the left and right graphs represent shapes corresponding to extreme negative ( −) and positive ( +) PC scores. Symmetric component (a) of PCA formed Tibet and WSP-HDM as distinct groups with some overlap, while the mix group was scattered between the two lineages. Along PC1, the change of leaf shape from subelliptical to suborbicular form was primarily related to the shape of the apical and basal regions and the length of petiole. The variation along PC2 mainly associated with the position of the maximum width of the leaf. For asymmetric component (b), all the specimens overlapped almost completely and the three lineages could not be discriminated, as has already been reported in several previous studies (Viscosi and Cardini 2011; Viscosi 2015; Liu et al. 2018). In detail, the variation along PC1 was mainly the changes in the relative position of the left/right sides at the maximum width of the leaf blade and the differences in the relative sizes of the leaf blade, whereas the variation along PC2 principally focused on the bending direction of the leaf blade toward left or right in asymmetric component

Fig. 7

Results of discriminant analysis (DA) of the shapes of leaves of the WSP-HDM vs. Tibet lineages (a), mixed individuals vs. Tibet lineage (b), and mixed individuals vs. WSP-HDM lineage (c). Black bars, Tibet lineage; red bars, WSP-HDM lineage; blue bars, mixed individuals