Geometric morphometric analyses of leaf shapes in two sympatric Chinese oaks: Quercus dentata Thunberg and Quercus aliena Blume (Fagaceae)

Abstract

Key message

Geometric morphometric analyses (GMMs) of the leaf shape can distinguish two congeneric oak species Quercus dentata Thunberg and Quercus aliena Blume in sympatric areas.

Contexts

High genetic and morphological variation in different Quercus species hinder efforts to distinguish them. In China, Q. dentata and Q. aliena are generally sympatrically distributed in warm temperate forests, and share some leaf morphological characteristics.

Aims

The aim of this study was to use the morphometric methods to discriminate these sympatric Chinese oaks preliminarily identified from molecular markers.

Methods

Three hundred sixty-seven trees of seven sympatric Q. dentata and Q. aliena populations were genetically assigned to one of the two species or hybrids using Bayesian clustering analysis based on nSSR. This grouping served as a priori classification of the trees. Shapes of 1835 leaves from the 367 trees were analyzed in terms of 13 characters (landmarks) by GMMs. Correlations between environmental and leaf morphology parameters were studied using linear regression analyses.

Results

The two species were efficiently discriminated by the leaf morphology analyses (96.9 and 95.9% of sampled Q. aliena trees and Q. dentata trees were correctly identified), while putative hybrids between the two species were found to be morphologically intermediate. Moreover, we demonstrated that the leaf morphological variations of Q. aliena, Q. dentata, and their putative hybrids are correlated with environmental factors, possibly because the variation of leaf morphology is part of the response to different habitats and environmental disturbances.

Conclusion

GMMs were able to correctly classify individuals from the two species preliminary identified as Q. dentata or Q. aliena by nSSR. The high degree of classification accuracy provided by this approach may be exploited to discriminate other problematic species and highlight its utility in plant ecology and evolution studies.

1 Introduction

There are many phenotypic traits which can be used for the differentiation of taxa, such as leaf morphology, anatomical descriptors, and other structural characters. Among them, leaf morphology is important for differentiating plant taxa (Costa et al. 2009), especially congeneric woody species, which tend to have little variation in floral features and short flowering seasons (Elias 1980). This is particularly true for oaks (Quercus L., Fagaceae), in which leaf features are “the most important” characters for discrimination (Stace 1982). Quercus is a genus of evergreen or deciduous shrubs and trees, and the largest genus in the Fagaceae family. About 450 species are widely distributed in temperate and subtropical regions of the northern hemisphere (Nixon 1993, 1997; Huang et al. 1999a, b; Manos et al. 1999), and 35 species are widely distributed in China (Huang et al. 1999a, b). Quercus taxonomy is extremely complex due to wide geographical distribution and species overlap, high levels of variability and phenotypic plasticity, and widespread hybridization where species ranges overlap (Rushton 1993; Bacilieri et al. 1996; Manos et al. 1999; Petit et al. 2005; Gerber et al. 2014).

Thus, discriminating oak species in mixed stands is not straightforward. However, progress has been made in the development of leaf morphological methods capable of resolving these problems in recent years. Traditional morphological measurements typically apply statistical techniques to a wide range of measurements, such as distances and distance ratios, angles, areas, and volumes (Mitteroecker and Gunz 2009). It can be efficient for discriminating the species. For example, Kremer et al. (2002) distinguished two oak species in nine mixed oak stands using 14 shape-related traits. However, traditional morphological methods are not effective for shape analysis, and do not allow interpretable graphic representation. Modern geometric morphometric methods (GMMs) using landmark and outline data are suitable for differentiating complex taxa because they can partition total coordinate variation into shape and non-shape (translation, rotation, and size) components, regardless of location, orientation, and scale factors (Slice 1996; MacLeod and Forey 2002; Adams et al. 2004; Mitteroecker and Gunz 2009; Viscosi et al. 2009). GMMs are also helpful for studying bilateral asymmetry and allometric patterns. Separate analysis of symmetric and asymmetric components of leaf variation in a principal component analysis (PCA) and two-block partial least-square (2B-PLS) analysis enables leaf shape variability to be recovered, thereby reducing the variation caused by developmental noise and instability (Viscosi 2015). Using GMMs, geometric information on morphological differences can be displayed in various images and diagrams, thereby highlighting differences in patterns and increasing the interpretability of results (Rohlf 2010). Due to these advantages, GMMs are increasingly used in morphological analysis (e.g., Viscosi and Cardini 2011; Viscosi and Fortini 2011; Souza et al. 2012; Viscosi et al. 2009; Tucić et al. 2018). Moreover, recent advances have shown that combinations of leaf morphological characters and molecular markers can significantly improve taxonomic delimitations (Viscosi et al. 2012; Rellstab et al. 2016).

In recent years, molecular markers were widely used in species discrimination due to the rapid development of DNA technology (e.g., Jaramillo-Correa et al. 2008; Du et al. 2011 for gymnosperm-conifers; Arnold et al. 2012; Eaton et al. 2015 for angiosperms). Various molecular markers have been used for species classification. Among them, nuclear microsatellite markers or simple sequence repeats (nSSRs) are bi-parentally inherited and experience extensive gene flow through both pollen and seeds (e.g., Petit et al. 2005; Du et al. 2013). Due to their co-dominant inheritance, high degree of polymorphism and relative ease of transfer between closely related species, nSSRs have become one of the most popular sources of molecular markers (Guichoux et al. 2011). Leaf shape is mainly genetically fixed but also affected by environmental factors (Graham et al. 2010). The environmental heterogeneity such as elevation, temperature, water, etc. can also affect the shape and size of leaves, e.g., leaves must be as wide as possible to absorb sufficient light and as flat and thin as possible to facilitate gas exchange, but will be dried out if too wide and thin (Tsukaia 2005). Thus, it is important to take environment factors into account when trying to explain the morphological differences.

Quercus dentata Thunberg and Quercus aliena Blume are generally sympatrically distributed in warm temperate forests, and they are important components of both deciduous broad-leaved and mixed broad-leaved forests (Ren and Wang 1985; Peng et al. 2007). The two species are found, and play major ecological roles, on southern slopes of mountains in northern China, and the neighboring countries Japan and Korea (Kazue and Atsuko 2015). In addition, like other oaks, the two species have great economic values, e.g., their leaves are used for sericulture, fruits for starch extraction, and bark for tannin extraction (Peng et al. 2007). Genetically, the two species are members of the Quercus (white oak) section, closely related, and regarded as sister species (Hubert et al. 2014). There are several major morphological differences between leaves of Q. dentata and Q. aliena: the former have shorter petioles, more densely stellate tomentose abaxial surfaces, more rounded bases, larger apices with shorter, blunter tips, and less secondary veins on each side of midvein compared with Q. aliena (Huang et al. 1999a, b). However, intermediate types of leaves have been observed in sympatric forests of Q. dentata and Q. aliena, which hinders discrimination of the two species (Huang et al. 1999a, b). Many previous studies on these two closely related species have mainly focused on their physiology and ecology (e.g., Fujiwara and Harada 2015; Box and Fujiwara 2015; Nugroho et al. 2016), while relatively few articles have addressed their morphological characteristics.

In this study, we collected numerous leaves from Q. dentata and Q. aliena trees and then investigated their morphological characteristics using selected leaf characters in a multivariate setting based on an a priori classification, to address the following three questions. First, can the two species be distinguished using selected leaf characters? Second, are symmetric and asymmetric patterns of variation helpful for distinguishing the two species and their hybrids? Finally, are leaf shape and size influenced by their environments, and if so which environmental factors are involved?

2 Materials and methods

2.1 Sampling strategy

We sampled trees in seven sympatric populations of Q. dentata and Q. aliena, covering most of their sympatric distribution in northern China (Fig. 1). Details of each population are listed in Table 1. Five intact, mature leaves were collected from each of 367 selected trees (which were spaced at least 5 m apart to minimize risks of selecting clones) for morphological analyses and DNA extraction. Thus, 1835 leaves were collected in total, all of which were taken to the laboratory for analysis, and voucher specimens were preserved at the Molecular Ecology Laboratory of Beijing Forestry University, China.

Fig. 1

Map showing the Quercus dentata and Quercus aliena sampling locations in this study, XS (XiangShulin), ZL (ZhongLou), DG (DongGou), DD (DaDianzi), QD (QiaoDuan), KY (KunYushan), and ZW (ZiWuyu)

Table 1 Locations and numbers of individuals of Quercus dentata and Quercus aliena based on genetic assignment

2.2 Genetic classification of the specimens

Genomic DNA was extracted from leaves of the 367 oak trees using a plant genomic DNA kit (Tiangen Biotech, China) following the manufacturer’s instructions. DNA concentrations in the extracts were measured by an ultramicro-spectrophotometer (Thermo Fisher, USA), and adjusted to a final concentration of 20–30 ng/μl for subsequent use. Twelve pairs of polymorphic nuclear simple sequence repeats (nSSRs) primers were selected for genotyping. These markers were originally developed in other oak species (Dow et al. 1995; Ueno et al. 2008, 2009; Durand et al. 2010) and the existence of repeats and polymorphisms were confirmed in the species we studied by Lyu et al. (2018). Detailed information about the primers, amplification, and genotyping procedures is presented in Supplementary Table S1. Bayesian cluster analysis was used to assign individuals to K clusters (1–10) on the basis of genotypes by STRUCTURE V2.3 (Pritchard et al., 2000). In each case, the program was set to run 200,000 Markov Chain Monte Carlo cycles (MCMC), following 100,000 burn-in cycles without any species identification information (USEPOPINFO = 0). Twenty iterations were performed for each K to estimate the best K according to Evanno et al. (2005). The most likely number of clusters was defined by computing ΔK and mean LnP(K) according Evanno et al. 2005 and Janes et al. 2017 (Fig. S1) using the structure harvester program (Earl and Vonholdt 2012). In order to validate the number of genetic clusters, we have provided two STRUCTURE plots (K = 2 and 3) for a visual comparison in Fig. S2. Admixture coefficient (Q) values were used to determine whether there were pure species or hybrids. The choice of the threshold of the Q value was dependent on efficiency and accuracy of assigning individuals to hybrid and purebred groups (Vähä and Primmer 2006). The reported threshold Q value for oaks is Q ≤ 0.9 (Lepais et al. 2009; Peñaloza-Ramírez et al. 2010) or Q ≤ 0.8 (Zeng et al. 2010). In this study, a strict threshold of 0.9 was chosen for a priori classification of the two species based on previous work in oaks (e.g., Lepais et al. 2009; Viscosi et al. 2012; Lyu et al. 2018). Individuals with Q ≥ 0.9 were considered to be purebreds, while individuals with Q < 0.9 from two genetic groups were considered to be hybrids individuals between two oak species.

2.3 Landmark configuration

The leaves were pressed, dried, and then scanned using a CanoScan 5600 F scanner (Canon Inc. Japan) at 600 dpi resolution. Scanned images were digitized to landmarks by Image J (Abràmoff et al. 2005). Thirteen landmarks were selected on each leaf, including the primary veins, symmetric structures, and homologous part (Fig. 2) (Jensen et al. 2002; Savriama and Klingenberg 2011; Viscosi et al. 2009). The landmarks consisted of three distributed along the middle axis of leaves (LM 1–LM 3) and nine others symmetrically distributed in other parts of the leaves (LM 4–LM 13). These landmarks were converted to a configuration of 13 pairs of Cartesian coordinates (x, y) for each leaf and the data were stored in “.txt” file format using the Image J program.

Fig. 2

Configuration of Quercus dentata and Quercus aliena leaves, showing locations of the 13 features used as landmarks in this study, with descriptions of the landmarks on the right

2.4 Statistical analysis

Cartesian x and y coordinates of the landmarks were used for the following analysis using the MorphoJ program (Klingenberg 2011). All the configurations were firstly rotated to generate a mean leaf configuration with minimal differences between corresponding landmarks by generalized Procrustes analysis (GPA) and were generated to a mirror image (Gower 1975; Rohlf and Slice 1990; Klingenberg et al. 2002; Klingenberg 2003; Zelditch et al. 2012). Secondly, leaf-level shape information was extracted for deviation detection by using the Procrustes distance to measure the absolute magnitude of the shape deviation and the Mahalanobis distance to provide an indication of the individual abnormality relative to others in the sample and a new dataset with no outliers (specimens that deviate strongly from the average) was constructed (Klingenberg and Monteirob 2005). Then, subsets of landmarks for the symmetric and asymmetric components of leaf shape were generated by GPA of the original configurations and their mirror image. Here, symmetrical component is composed of median axis and joint mirror image shifts of pairs of corresponding landmarks in any direction, while asymmetric component quantifying differences between the original configurations and mirrored images as described in the manual of MorphoJ (Klingenberg et al. 2002). Fourthly, the variation of leaf shape hierarchically distributed among leaves, trees, and species was used to create mean configurations at each of these levels. Leaf-level data were used for leaf shape variation analysis, tree-level data for species discrimination, and species-level data (including the hybrid individuals) for multivariate analysis. Centroid size, symmetric, and asymmetry components of leaves were used for further analysis. Finally, matrix correlations of each hierarchical level were used as datasets for subsequent multivariate statistical analyses.

The normalized matrices were then subjected to PCA to identify latent (principal) components capturing most of the variation in the dataset (Benzécri 1992), i.e., global variation in leaf landmark configurations, and identify the leaf shape variability at leaf level. To investigate the leaf morphology variations of the entire samples, the symmetric and asymmetric components were separately subjected to leaf-level PCA. Because leaf shape variability (symmetric and asymmetric components) might be associated with leaf size (i.e., allometry) (Klingenberg 1998), 2B-PLS analysis was performed to detect covariation of leaf shape and size with a 10,000 permutation test using MorphoJ.

Discriminant analysis (DA) and canonical variate analysis (CVA) were used for tree-level species discrimination. Generally, DA is the most widely used statistical method for investigating taxonomic differences, focusing on two groups, while CVA is used for investigating three or more groups (Klingenberg 2011). In this study, DA was performed using cross-validated scores classification tables with T² statistics (P value for tests with 1000 permutations < 0.0001) while CVA was performed with permutation tests using Mahalanobis distances (T² statistics; 10,000 permutations per test). Results of the DA (which compares each pair of groups separately) were visualized in a frequency bar chart and the CVA results were presented in a scatter plot of the significant canonical variate (CV) scores.

GMMs allow shape differences between specimens to be recognized, using diagrams to visualize variability along the axes (i.e., PC/CV/PLS). In this study, an overlapping of transformation grid and wireframe graph was used to reconstruct the leaf shape by thin-plate spline analysis (Bookstein 1996; Klingenberg 2011; Viscosi and Cardini 2011; Viscosi 2015).

Linear regression analysis was used to study the relationships between environmental variables and both leaf shape (symmetric and asymmetric components) and size parameters. The leaf shape and size dataset was generated by MorphoJ, using symmetric and asymmetric components that contributed most strongly to PC1 and thus explained most of the variation in leaf shape (61.7 and 61.6%, respectively). The climate data were generated using bioclimatic parameters from the WORLDCLIM database (http://www.worldclim.org) for the current climate across the collection localities at 30 s resolution by ArcMAP 10.2 (http://www.esri.com/software/arcgis/arcgis-for-desktop). Because collinearity can lead to instability of parameter estimates in regression models (Chatterjee et al. 2006; Dormann et al. 2013), environmental factors that were strongly correlated with other environmental factors were removed by variance inflation factor (VIF) analysis with the threshold values of 0.7 implemented in the “usdm” package, R. Finally, mean temperature of the warmest quarter (MTWQ), precipitation in the wettest month (PWM), and altitude were selected from 19 climate factors and three geographical factors (longitude, latitude, and altitude) for the analysis of correlations between leaf parameters and environmental variables. Values of the environmental factors at each population’s location are listed in Table S2.

3 Results

3.1 A priori assignment of the individuals to the two species

STRUCTURE analysis clearly divided the 367 individuals into two main clusters with different genotypes, with indications of wide and asymmetric gene flow between them (Fig. S1). The genetic analyses assigned 129 individuals to pure Q. dentata, 194 individuals to pure Q. aliena, and 44 individuals to hybrids, with a 0.9 admixture coefficient (Q) threshold (Table 1).

3.2 Leaf morphology analysis

3.2.1 Multivariate analysis of raw data

The raw coordinate matrix obtained from image analysis of the 1835 leaves was translated into a normalized matrix by Procrustes fitting, and the configurations were rotated to the concentrated distribution around 13 leaf landmarks based on the raw coordinate matrix (Fig. 3a), symmetric components (Fig. 3b), and asymmetric components (Fig. 3c). After removing outliers from the landmark data matrices, data for 1815 of the 1835 leaves were retained and used for subsequent multivariate analyses.

Fig. 3

Results of generalized Procrustes analysis of the leaf shape of Quercus dentata, Quercus aliena, and hybrid based on the full raw coordinate matrix (a), symmetric components (b), and asymmetric components (c)

3.2.2 Interspecific leaf shape variation

In the scatter plot of PC1 and PC2 scores for symmetric components generated by PCA, leaves of the two species formed distinct groups with some overlap, whereas leaves of their hybrids were scattered (upper panel of Fig. 4a). In detail, the first and second PCs accounted for more than two thirds of the total leaf shape variation (61.7 and 14.6%, respectively) and the first five of PCs accounted for 96.6% of the total leaf shape variation (Table S3). Symmetric components with extreme PC1 and PC2 scores (the maximum and minimum values of the scatters projected on both ends of the PC axis) revealed the main differences between leaves of the two species (lower panel of Fig. 4a). The PC1 results showed that Q. aliena leaves had a longer petiole (distance between LM 1 and 2); narrower blade tip (manifested by contraction of the grid and squeezing between the leaf tip at LM 3 and widest point, at LM 6 and 11, relative to the maximal width), and relatively wider basal region (distance between LM 8 and 13) compared to Q. dentata leaves. The distance between LM 7 and 12 was related to the degree of leaf blade lobes division. PC2 was mainly related to the locations of lobes at the maximal width (LM 11 and 6) and their sinuses (LM 12 and 7), the negative values were associated with deep lobes, and positive values were associated with shallow lobes.

Fig. 4

Results of leaf-level principal component analysis (PCA) of Quercus dentata, Quercus aliena, and hybrid leaves’ symmetric components (a) and asymmetric components (b). Scatter plots of PC1 and PC2 scores, with 95% confidence ellipses in a. Transformation Grids: the left and right graphs represent shapes corresponding to extreme positive (+) and negative (−) PC scores

In the PC1 versus PC2 scatter plot for asymmetric components, the specimens were rather fully and densely overlapping than scattered (Fig. 4b). In detail, the first and second PCs accounted for 76.1% of the total leaf shape variation (61.6 and 14.5%, respectively): PC1 reflected the asymmetry of both lobes at the greatest width of the leaf while PC2 reflected the asymmetry due to bending of the leaf blade.

3.3 Allometry

2B-PLS analyses of relations between allometry and symmetric components at the tree level resulted in a RV coefficient (squared Pearson correlation coefficient) of 0.1089, P < 0.001 (Fig. 5a), indicating significant covariation between these components and leaf size. Graphical reconstruction of the symmetric components along PLS2 showed that negative values are associated with longer petioles and narrow to obovate leaf shape, while positive values are associated with short petioles and narrow to inverted water drop leaf shape. The graphical reconstruction along PLS1 (log centroid size axis) indicated that Q. dentata has larger leaves and positive leaf shape values while Q. aliena has smaller leaves with negative leaf shape values. However, no significant allometric patterns of asymmetric variation were detected for leaves of Q. aliena, Q. dentate, or their hybrids (RV coefficient, 0.0068, P < 0.1673) (Fig. 5b).

Fig. 5

Scatterplots obtained from 2B-PLS analysis of the relationship between size (log centroid size) and both symmetric components (a) and asymmetric components (b) of the shape of leaves of Quercus dentata, Quercus aliena, and their hybrids. 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 components

3.4 Species discrimination

The cross-validated DA revealed significant differences in leaf shape between Q. dentata vs. Q. aliena (Fig. 6a), large overlap in shapes of leaves of Q. aliena vs. the hybrids (Fig. 6b), and small overlap in shapes of leaves of Q. dentata vs. the hybrids (Fig. 6c). For the pairwise comparisons, T² values (all significant, P < 0.0001) were 1659, 220, and 84 for Q. aliena vs. Q. dentata, hybrids vs. Q. dentata and hybrids vs. Q. aliena, respectively. In these comparisons, 95.9 and 96.9% of Q. aliena and Q. dentata, 64.4 and 81.5% of the hybrids and Q. aliena, and 71.1 and 97.7% of the hybrids and Q. dentata trees were correctly classified, respectively.

Fig. 6

Results of discriminant analysis (DA) of the shapes of leaves of Quercus dentata vs. Quercus aliena (a), Q. aliena vs. hybrids (b), and Q. aliena vs. hybrids (c). Black bars, Q. dentata; red bars, Q. aliena; blue bars, hybrids

Scatter plots of CVA scores showed that leaves of Q. dentata and Q. aliena trees mainly distributed along the CV1 axis (61.7%) with scores of − 4 to 0 and − 2 to 4, respectively (Fig. 7). In addition, the leaf morphology of the hybrids is intermediate to some extent but also highly variable with a certain overlap with Q. aliena. Overlapping transformation grid and wireframe graphs showed that the main differences between leaves of Q. dentata and Q. aliena lay in the patterns of expansion and contraction from the leaves’ bases to their apexes (Fig. 7). In contrast, Q. dentata, Q. aliena, and hybrid leaves largely overlapped along CV2 (38.2%).

Fig. 7

Results of tree-level canonical variate analysis (CVA) of Quercus dentata, Quercus aliena, and hybrids’ leaf shapes (symmetric components) showing 90% confidence ellipses for the species. Transformation grids represent leaf shapes reflecting the negative (− shape) and positive (+ shape) extremes of the CV1 axis

3.5 Regression analysis of relations between environmental variables and leaf shape and size parameters

The linear regression analyses showed clear correlations between the leaf morphology and environmental factors (Table 2). Symmetric components of Q. aliena leaves were significantly correlated with MTWQ, PWM, and altitude, those of Q. dentata leaves were correlated with MTWQ and altitude, while those of hybrids were significantly correlated with PWM. In contrast, asymmetric components of both sympatric oak species and their hybrids were only correlated with MTWQ. In addition, leaf size of Q. aliena was significantly correlated with altitude, leaf size of Q. dentata was significantly correlated with PWM, and leaf size of the hybrids was correlated with MTWQ and PWM.

Table 2 Results of linear regression analysis of relationships between environmental factors

4 Discussion

GMMs based on a priori classification of genetic analysis can provide reliable approaches for quantitative and graphical resolution of variation in leaf shape, with several major advantages over traditional morphological analysis (Viscosi et al. 2012). Based on the a priori classification by nSSR and GMMs, we were able to develop a robust methodology for discriminating the sympatric oak species and explore the relationships between leaf allometry, shape, size, and environmental factors, which might be useful for the morphological studies of closely related oaks in China.

The results showed that analysis of symmetric components of leaf shape can successfully distinguish the two species if we presume that nSSR assignment well reflects taxonomic separation. In PCA score plots for symmetric components, Q. dentata and Q. aliena leaves formed two clearly separated groups, while leaves of their hybrids were scattered (Fig. 4a). Deformation grids showed that Q. dentata leaves have a shorter petiole, narrower basal region, deeper clefts, larger apical region, less lobed blade, and broader blade tip than Q. aliena leaves. Moreover, the morphological discrimination analysis after a priori classification at the molecular level proved to be useful to distinguish the two species and their hybrids. The DA analysis highlighted significant differences in leaf morphology variables between Q. dentata and Q. aliena (Fig. 6). Leaf shape DA could correctly classify genetically pure individuals with > 95% probability, but only correctly classify hybrids with much lower probability (ca. 64–71%). Interestingly, DA provided higher probability to correctly distinguish Q. dentata from hybrids than Q. aliena from hybrids. In addition, in CVA score plots, the two species were significantly separated along CV1, and the hybrids were mainly located between them, albeit closer to the Q. aliena cluster (Fig. 7). The CVA and DA results suggest that the hybrids’ leaf shapes were affected by both parental species, but most strongly by the dominant species, Q. aliena. Interestingly, a corresponding pattern of variation in hybrids and parental species has been observed in populations of European white oaks (Quercus frainetto and Quercus petraea) (Viscosi et al. 2012; Stephan et al. 2018). This phenomenon may be caused by extensive and imbalanced gene flow between oak species (Gerber et al. 2014) but needs further investigation.

Several studies have found that symmetric components of leaf shape are related to genetic factors, while asymmetric components are generally related to environmental perturbations (Yoshioka et al. 2004; Bresson et al. 2011; Viscosi 2015). However, both the symmetric and asymmetric components of the shape of leaves studied here were significantly related to environmental factors. Moreover, responses of the symmetric components of leaves of the two species and their hybrids to several environmental factors differed, while the only tested environmental factor that influenced the asymmetric components was the mean temperature of the warmest quarter, MTWQ. It is reported that leaf phenology and ecophysiological properties are more environmentally influenced than genetically controlled in highly adaptive and plastic genotypes as described in Vitasse et al. (2010) and Bresson et al. (2011). However, further investigations such as controlled experiments are needed for the validation of leaf shape adaptation to different environments.

The size of Q. aliena leaves was significantly related to elevation (Table 2) as reported for other oaks (Viscosi 2015; Bresson et al. 2011). However, the leaf sizes of Q. dentata and hybrids were related to precipitation in the wettest month or mean temperature of the warmest quarter. Thus, leaf size can clearly be influenced by different environmental factors and genotypes, and their effects warrant further attention.

In conclusion, the GMM analysis of the leaf morphology of the two sympatrically distributed oak species suggests that we can clearly separate leaves of Q. dentata and Q. aliena morphologically. Moreover, individuals classified as putative hybrids through genetic analysis were consistent with morphologically intermediate individuals. In field observations of Q. dentata and Q. aliena at the XS and ZL sites in 2017, we found that their flowering times were asynchronous but the male florescences of Q. dentata and female florescences of Q. aliena overlap phenologically to a limited extent. The overlapping of the flowering time of the species suggests the possibility of introgression or hybridization between the two species. Introgression as an adaptive process is a common phenomenon in some tree species, such as oak and poplar (e.g., Lepais et al. 2009; Suarez - Gonzalez et al. 2018). But carefully controlled pollination experiments are required to confirm this possibility especially on the non-sympatric areas of the two species.

5 Conclusion

Thorough understanding of variation in leaf shape is highly important for taxonomists and ecologists, and morphological variation of Quercus components in mixed forest stands has interested biologists and botanists for decades (Kremer et al. 2002; Gugerli et al. 2007). GMMs based on genetic classification provided detailed, accurate information on the morphological variation in two sympatrically distributed Chinese oaks, and revealed correlations between leaf morphology and environmental factors. The data we have acquired could facilitate research into the taxonomy, phylogeny, and ecology of oaks. Moreover, extending the applied approach to studies of other plants will likely provide new insights into the biology and evolution of both oaks and other taxa.