Morphological variation in mixed oak stands (Quercus robur and Quercus petraea) is stable across western European populations

Antoine Kremer (1), Jean Luc Dupouey (2) , J. Douglas Deans (3),  Joan Cottrell (4), Ulrike Csaikl (5), Reiner Finkeldey (6), Santiago Espinel (7), Jan Jensen (8), Jochen Kleinschmit  (9), Barbara Van Dam (10), Alexis Ducousso (1), Ian Forrest (4), U. Lopez de Heredia (11),   Andrew J. Lowe (3), Marcela Tutkova (5),  Robert C. Munro (3),  Sabine Steinhoff (9), Vincent Badeau (2)

1 Institut National de la Recherche Agronomique (INRA), Unité de Recherches Forestières, BP 45, 33611 Gazinet cedex, France.
2 Institut National de la Recherche Agronomique (INRA), Unité d’Ecophysiologie Forestière, 54280 Champenoux, France.
3 Centre for Ecology and Hydrology, Edinburgh, (CEH), Bush Estate, Penicuik, Midlothian, EH26 0QB, Scotland, U.K.
4 Forestry Commission (FC), Forest Research, Northern Research Station, Roslin, Midlothian, Scotland, U.K., EH25 9SY.
5 Austrian Research Centre (ARCS), Seibersdorf. A-2444, Austria
6 Swiss Federal Research Institute WSL, Zürcherstrasse 111, CH-Birmensdorf, Switzerland.
7 NEIKER A.B., Granja Modelo-Arkaute, 01080 Vitoria-Gasteiz, Spain
8 Danish Forest and Landscape Research Institute (DFLRI), Hørsholm Kongevej, Hørsholm DK-2970, Denmark.
9 Niedersächische Forstliche Versuschsanstalt Fortamstr.6, D-34355 Staufenberg-Escherode, Germany.
10 ALTERRA, Green World Research Institute, Wageningen, The Netherlands
11 Present address: ETSI Agrarias, Universidad de Valladolid, Avda Madrid 57, 34004 Palencia, Spain

Abstract
Leaf morphology was assessed in nine mixed oak stands (Quercus petraea and Q. robur) located in eight European countries. Exhaustive sampling was used in an area of  each stand where the two species coexisted in approximately equal proportions (about 170 trees/species/stand). Fourteen leaf characters were assessed on each of 5 to10 leaves collected  from the upper part of each tree. Three multivariate statistical techniques (CDA, Canonical Discriminant Analysis; PCA, Principal Component Analysis; MCA, Multiple Correspondence Analysis) were used in two different ways: first on the total set of leaves over all stands (global analysis) and second, separately within each stand (local analysis). There was a general agreement of the results among the statistical methods used and between the analyses conducted (global and local). The first synthetic variable derived by each multivariate analysis exhibited a clear and sharp bimodal distribution, with overlapping in the central part. The two modes were interpreted as the two species, and the overlapping region was interpreted as an area where the within-species variations were superimposed. There was no discontinuity in the distribution or no visible evidence of a third mode which would have indicated the existence of a third population composed of trees with intermediate morphologies. Based on petiole length and number of intercalary veins, an “easy to use “ discriminant function applicable to a major part of the natural distribution of the species was constructed. Validation on an independent set of trees provided a 98% rate of correct identification. The results were interpreted in the light of earlier reports about extensive hybridization occurring in mixed oak stands. Maternal effects on morphological characters, as well as a lower frequency or fitness of hybrids in comparison with parent species could explain the maintenance of two modes , which might be composed of either pure species or pure species and introgressed forms.
Keywords: morphology, leaf, Quercus robur, Quercus petraea, taxonomy, Europe, mixed stand

1 - Introduction
	Morphological variation in mixed oak stands composed of Quercus petraea and Quercus robur has been of general interest in Europe for decades. The two species commonly coexist  in mixed stands and foresters need rapid “easy to use" morphological screening methods which would  help to discriminate between the species (Cousens, 1963; Carlisle and Brown, 1965; Granjean and Sigaud, 1987). Dendrologists require classification criteria for taxonomic purposes (Olsson, 1975a; Kissling, 1980a and b; Minihan and Rushton, 1984; Aas and Friedrich, 1991). Forest biologists and ecologists  seek unambiguous traits  which could be used for studying introgression between these two interfertile species (Wigston, 1975; Rushton, 1978, 1979; Iestwaart and Feij, 1989). The vast amount of literature devoted to morphological variation in mixed oak stands demonstrates the debate that has been raised. In addition to efforts  to  define the most appropriate morphological traits to be assessed (Olsson, 1975b, Kissling, 1977; Rushton, 1983; Safou and Saint Martin, 1989; Aas, 1995), investigations  have been conducted on  relevant statistical methods (Dupouey, 1983, Dupouey and Badeau, 1993; Kleinschmit et al., 1995). However these reports have not led to  general agreement on the methods and traits to be used for distinguishing Quercus petraea from Q. robur, mainly because no comparative studies were conducted across different countries. Key issues in the controversy were also the sampling and analytical methods for assessing morphological variation between two interfertile species, and particularly, the inclusion or not, of trees exhibiting intermediate "morphologies". In several cases, investigations were made on trees that were a priori selected as belonging to groups with contrasting morphologies, purposely excluding ambiguous trees. Sampling procedures or data analyses based on a priori grouping inevitably introduced circularity in the taxonomic assignment procedure. Alternative methods for species identification have been explored in recent years using either isozymes (Zanetto et al., 1994; Bacilieri et al., 1995; Samuel et al., 1995) or molecular techniques (Moreau et al., 1994; Bodénès et al., 1997; Muir et al., 2000). These methods generally assume that species separation was known and were used to compare frequencies of genetic markers between the two species. In these studies, the original information used for species separation was leaf or fruit morphology, and we were back to the initial problem. Even for molecular markers to be used, there is a need for clarification of species distinction between Q. petraea and Q. robur based on morphological traits. 
Our contribution is an attempt to provide a unifying method to assess morphological variation for species recognition that would be applicable across countries and forests. This investigation  formed a necessary step in a project aimed at comparing gene diversity between the two species in different European countries. As there was no general agreement on a method for taxonomic assignment of temperate oaks, we combined protocols available in the literature in order to design a standard procedure  which could be  applied in widely separated forests for taxonomic assignment  prior to the gene diversity investigation. With this general aim, we proceeded in three successive steps. First we compared multivariate statistical methods for the discrimination of trees according to their leaf morphology. At this stage, there was no a priori classification of the oak trees in two or three groups. Only trees were the units to be separated. In a second step, we verified  that the morphological variation of trees based on the statistical methods used was independent of geographical region. The procedure was therefore based on samples originating from widely separated stands throughout western Europe. Lastly we constructed a discriminant function that could be applied as an operational “easy to use” procedure over a  wide geographic range for species recognition. Finally, this function was tested on a totally independent set of trees.
2 – Material and methods
Study stands 
Nine mixed stands  of Q. petraea and Q. robur were sampled in 9 regions of Western Europe. The criteria used in the selection of the stands were age, composition and size. Only adult stands  of several hundred trees were selected. The stands were chosen so that the two species were present in approximately equal proportions, if possible.  All adult trees within a delineated area of the stand were used for the morphological assessment: the sampling was exhaustive, with no a priori selection of trees. The exception to this rule was the stand of Salinasco, where local constraints did not permit studies to proceed as an exhaustive sampling.  All other stands contained on average  369 oak trees within the study  area (Table 1). In total 3025 trees were used for the leaf morphological assessment.
Assessment of leaf morphological traits
 Up to 10 (usually 5) fully expanded leaves were sampled in the mid- to upper crown of the trees. As far as was practicable, they were insect and disease free,  and were collected  from the first flush of the year. In total 16055 leaves were measured. The protocol of leaf morphology assessment was based upon Rushton (1983) and Kissling (1977) with some modifications aimed at simplifying the procedure. The following variables were assessed on each leaf:
•	Five dimensional characters (Figure 1a) 
Lamina length (LL), petiole length (PL), lobe width (LW), sinus width (SW), length of lamina at largest width (WP). WP and LW were measured at the tip of the widest lobe of the leaf. 
•	Two counted variables 
Number of lobes (NL): total number of lobes including those on the right and the left part of the leaf, except the terminal lobe (irrigated by the midrib). A lobe was considered to be  present when it was clearly irrigated by an axillary vein. There was no size limit to a lobe.
Number of intercalary veins (NV): an intercalary vein was a secondary vein irrigating a sinus and extending at least half way from the midrib to the base of the sinus.
•	Two observed variables 
Basal shape of the lamina (BS): this was scored as an index varying from 1 to 9 as indicated on Figure 1b. 
Abaxial laminar pubescence (PU): a scoring of the density of pubescence  was made using Kissling’s grading system from 1 (no pubescence) to 6 (dense hairiness) (Kissling, 1977). Pubescence was assessed with a stereomicroscope (x 30) and took into account both stellate and simple hairs. 
•	Five transformed variables
Lamina shape or obversity (OB): 	OB = 100* WP / LL
Petiole ratio (PR): 			PR = 100 * PL / (LL+PL)
Lobe depth ratio (LDR): 			LDR = 100 * (LW-SW) / LW
Percentage  venation (PV): 		PV = 100 * NV / NL
Lobe width ratio (LWR): 		LWR = 100 * LW / LL 
Data analysis
Three different multivariate procedures were used for analyzing the data: (i) canonical discriminant analysis (CDA) of the 14 original variables, using the tree as the classification variable; (ii) principal component analysis (PCA) of the 14 variables; (iii) multiple correspondence analysis (MCA) of the 14 variables, each variable being divided into 15 classes of nearly equal weight (Benzécri, 1992). The three methods aimed at combining the original variables into independent synthetic variables that explained the greatest part of the total variation observed among the trees. The units that were discriminated were the trees: there was no a priori grouping into two or any other number of groups. The differences among the different methods were related to the criteria used to construct the synthetic variables:
•	In CDA, the synthetic variables (canonical variates) were linear combinations of the original variables, constructed so that the ratio (variance between trees / variance within trees (between leaves)) was maximal. Variables taking integer values only, such as counted variables, could introduce bias in the calculation because within tree variance was frequently null.  This was the case for pubescence, which displayed a null variance among leaves in half of the sample of trees. Thus, this variable was discarded from all CDA analyses.
•	In PCA, the synthetic variables (principal components) were also linear combinations of the original variables (mean values of the five leaves per tree) that were constructed so that they displayed the largest variance between trees. Computations were performed on the correlation matrix between original variables.
•	In MCA, the synthetic variables (principal axes of inertia or factorial axes (Benzécri, 1992)) also had the largest variance between trees. Computations were performed on the indicator matrix with trees as rows and categories of variables as columns. In contrast to PCA, the synthetic variables included non linear relationships between the original variables.
The three multivariate methods were then applied in two different ways:
•	on the overall data set using all 3025 trees of the 9 sites (global analysis). In this case, CDA was based on the 14 variables*16055 leaves, PCA on the 14*14 correlation matrix computed from the 14 * 3025 trees table, and MCA on the 14 variables * 3025 trees table, each variable being divided into 15 classes of nearly equal weight. 
•	on each separate study stand  (local analysis). The three statistical methods were applied separately within each stand by using only the data of that stand (but maintaining the same 15 classes boundaries for all variables in MCA across all sites).
3 - Results
Leaf size effect
Because the data originated from widely separated stands, it was suspected that the original variables might be affected by the size of the leaves. Hence the correlation between LL (Lamina length) and all other traits was computed at the leaf level over the whole data set (Table 2). As expected, the dimensional traits showed positive correlations, whereas all other traits were not influenced by size effects. There was however, an important variation between dimensional traits: LW and WP were strongly correlated with LL, whereas SW and PL were only moderately correlated. The transformation of the original variables, which consisted mostly in taking ratios between dimension traits, reduced the correlations as can be seen by the low correlations between OB, PR, LDR and LWR with LL. An interesting observation was that the morphological traits that were traditionally used in the literature for species separation (PL or PR, NV or PV, BS) showed only weak correlations with LL.
Multivariate analyses over the whole data set (global analysis)
The three methods were applied over all of the 3025 trees and the distributions of the synthetic variables were plotted for each method (Figure 2). In each method, the first synthetic variable contributed to a major part of the total variation (Table 3) as expected, but in addition,  the nextsynthetic variables 2nd and 3rd)  exhibited a much lower contribution. Whatever the statistical method used, the distribution of the first synthetic variable showed a  bimodal distribution as shown by Figure 2. The whole set of trees was therefore composed of two populations that were partially overlapping (Figure 2). There was a strong similarity  between the position of the trees  derived from the three methods, as shown by the rank correlation of their values for the first synthetic variable (Table 4). CDA provided slightly less congruent results.
Multivariate analyses within each stand (local analysis)
The same analyses that were conducted over the whole data set were applied to each individual stand. Again in each stand, whatever the statistical method used, the first synthetic variable showed a bimodal distribution (data not shown). The bimodality was less  apparent in Dalkeith and Hald Ege, mainly due to the unbalanced representation of the two species in these two stands, as shown later (Table 8).  To examine whether the distribution of the trees along the first synthetic variable in the local analysis was congruent with the distribution observed in the global analysis, the rank correlations between the first synthetic variables between the two analyses were computed (Table 5). These correlations remained extremely high, except for Hald Ege, suggesting that the analyses  for each stand provided the same distribution along the first synthetic variables as  the analyses conducted over the whole data set. 
Biological significance of the synthetic variables
The two analyses (over the whole data set, and within each separate stand ) and the three statistical methods, indicated that the sample of trees that was analysed in this study was composed of two populations that corresponded to the two modes observed  in Figure 2. There was neither discontinuity in the distribution nor an additional mode that would correspond to the existence of a third population. However, the two bell shaped distributions overlapped at their tails. It was more likely that the overlapping region encompassed the within population variation rather than corresponding to a third separate population. In order to  characterise the two populations, we focused on the biological significance of the first synthetic variables that discriminated the groups, by computing the correlations between the synthetic and original variables for the global analyses (Table 6). The original variables that exhibited the highest correlations with the first synthetic variable were the same for the three statistical methods used: petiole length (PL and PR), intercalary  venation (NV and PV), pubescence (PU) and sinus width (SW and LDR). The basal shape of the lamina (BS) and number of lobes (NL) contributed moderately to the first synthetic variable. Other morphological variables were only weakly correlated with the first synthetic variable. Interestingly, the original variables showing the highest correlation with the first synthetic variables were those that are traditionally used for species identification in oaks. Thus, the first synthetic variable could be interpreted as a gradient between Q. robur and Q. petraea There were some slight discrepancies between the three statistical methods in regard to the contribution of the original variables to the synthetic variables. MCA and PCA were fully congruent, but the first synthetic variable in CDA was characterized by a lower contribution of all variables except for PL and PR; the contribution of sinus width (SW or LDR) was especially low in comparison  with the other two methods.
The next synthetic variables (2nd and 3rd) showed a continuous distribution in contrast to the first one (data not shown). The second synthetic variables of MCA and PCA were mostly correlated with dimensional traits (LL, LW and WP) and were therefore interpreted as a leaf size gradient. In the case of CDA, the second synthetic variable was still correlated with traits that exhibited species differences. In summary, among the 14 original morphological variables, 7 were strongly correlated with the first synthetic variable (whatever the statistical method used), 3 were correlated with the second axis that corresponded to leaf size, and the remaining 4 were distributed among the remaining 12 synthetic variables. 
Operational method for species assignment
Because the three multivariate statistical methods provided the same bimodal distribution along the species gradient (Figure 2), and because these results were congruent across stands (Table 5), we attempted to construct an operational species assignment procedure, that had the following criteria (i) easy to use, (ii) applicable across a large geographic range, (iii) with  minimal misclassification.
The following "consensus" clustering procedure was used in order to define typical Q. robur and Q. petraea trees: for each of the three multivariate methods separately, trees were classified as either Q. robur or Q. petraea according to their position along the first synthetic axis, using the standard k-means clustering algorithm (procedure FASTCLUS of SAS, 1989). Then, trees classified in the same species by all three methods were retained as typical of this species, the other ones being considered as unclassified trees ( 4.7% of the total, Table 7). Once again, there was strong agreement between the three methods, CDA being slightly  distinct from MCA and PCA (Table 7). A stepwise discriminant analysis was performed using these two groups of "typical" trees as the classification variable. The first discriminant variables selected were PR, PV, PU and BS, with partial R2 of 64, 44, 13 and 16% respectively.  Only the first two “easy to assess” variables PR and PV were retained, but we substituted NV for PV and PL for PR (easier to measure and only slightly less efficient). Based on NV and PL, the following discriminant function for oak identification was proposed:
ID = 357 - (97 x PL) + (205 x NV)                   (PL in mm)
This function gave positive ID values for Q. robur and negative values for Q. petraea. The percentage of misclassifications was as low as 1% in the calibration set of "typical" trees (Table 7). Figure 3 shows the distribution of ID values among all the 3025 trees measured, and Figure 4 presents the position of Q. robur and Q. petraea trees (classified according to the consensus clustering procedure) as a function of NV and  PL values.
We used a resampling procedure to assess the  effect of the number of leaves on the ID values. Subsamples were built from the initial total sample by randomly reducing the number of leaves per tree to four, three, two or one. New ID values for each tree were computed using this number of leaves, and compared with ID values based on 5 leaves per tree. The average tree misclassification percentages over 1000 iterations of the random selection process increased slowly when the number of leaves taken into account decreased: 1.1%, 1.8%, 2.9% and 5.6% for four, three, two and one leaves per tree, respectively. However, the intra-specific order of ID values changed faster when the number of leaves decreased: whatever the species, the Spearman rank correlation with ID values based on 5 leaves per tree were 0.98, 0.95, 0.90 and 0.80, respectively.  A sample of three leaves per tree seemed to be an optimum compromise.
Finally, the discriminant function was validated using an  independent sample of 773 oaks from northeastern France (Dupouey & Badeau, 1993), a region not covered by the present study. In that study, trees had been classified as either Quercus robur, Quercus petraea, Quercus pubescens or intermediate based on a MCA analysis. Using our discriminant function and 5 leaves per tree, only 1.6% of the Q. robur and Q. petraea trees were not correctly classified. Interestingly, 98% of the Q. pubescens trees were classified as Q. petraea and intermediate trees were distributed between Q. robur (27%) and Q. petraea (73%).  
The proportion of trees belonging to the Q. robur, Q. petraea or unclassified groups (according to the consensus clustering procedure) varied according to the geographic location of the stands (Table 8). Whereas stands from Southern, Western and Central Europe displayed a very low proportion of unclassified trees (2% or less), those from The Netherlands, Denmark and Great Britain had between 5 and 13% unclassified trees. This higher proportion of unclassified trees in northern stands was due to a higher morphological similarity between trees from the two species in these stands than elsewhere in Europe. Q. petraea trees from northern stands had shorter petioles than in other parts of Europe (11.3 and 15.2 mm respectively), while Q. robur displayed  greater pubescence (3.6 and 1.4, respectively).

4. Discussion

	The three multivariate statistical methods clearly demonstrated that the 9 mixed oak stands that were analyzed in this study  consisted of 2 populations, as indicated by the bimodal distribution of the synthetic variables. A closer examination of the morphological characters included in the synthetic variables showed that the two modes in the double bell shaped curves corresponded to Q. petraea and Q. robur. The conclusions of this study were reinforced by the congruence of the results across sites.  Whether the data were analyzed separately within each site or on the total set of leaves over all sites, the same results were obtained: the bimodality was observed within each analysis and the distribution of the trees for the first synthetic variables were highly correlated among the local and global analysis (Table 5). The only exceptions to this were Hald Ege (in Denmark) and, to a lesser extent, Dalkeith (in Scotland), because of an unbalanced distribution of trees between the two species (in each case, more than 90% of the trees belonged to one species). The robustness of the results was also due to the congruence of the conclusions drawn from three separate multivariate  analyses (Table 4). There were only a few discrepancies among the three methods, particularly between CDA and the two other methods (MCA and PCA). Because CDA also included the within tree variation, it was less affected by the “noise” due to sampling of leaves within trees. As a result, among the three statistical methods tested, CDA provided slightly more congruent results between the global and local analyses (Table 5). These results were comparable with other investigations  using multivariate  analysis conducted in various parts of western Europe (Grandjean and Sigaud, 1987, in the Central part of France ; Dupouey and Badeau, 1993 in the Eastern part of France ; Iestwaart and Feij (1989) in the Netherlands), central Europe (Aas, 1993, in Germany and Poland) and eastern Europe (Borovics, 1999), although they differed in their sampling strategy. All these case studies also showed a strong bimodal distribution of the synthetic variables.
The sharp bimodality does not provide evidence for  the existence of a third population composed of intermediate phenotypes. The general picture of the total diversity of leaf morphology in the Q. petraea-Q.robur oak complex was rather the coexistence of two populations with some overlapping in their morphological distribution (Figure 2), rather than the coexistence of three populations, the additional one being composed of intermediate phenotypes. If a third population had to be defined it  would have to rely on some “ a priori ” limits that would be difficult to  identify  from the distribution  on the synthetic axis (Figure 1 ).  Alternatively, it could be based on morphological variables other than the ones we used in this study, although we took into account most of the leaf morphological variables known to discriminate between the two species. Our data rather suggested that trees exhibiting intermediate phenotypes corresponded actually to either Q. petraea or Q. robur. This raised the  question of species assignment in the overlapping region of the distributions (Figure 2). There was a region of uncertainty in the tails of the distributions, where species assignment was subject to error. A statistical procedure to estimate the error rate could be calculated after fitting theoretical distribution curves to the empirical distribution observed in Figure 2. An alternative way would have been to include additional diagnostic characters available for taxonomic identification. In a number of different reports intermediate phenotypes have received special attention, first by delineating a specific class of intermediate phenotypes in comparison to reference populations (Elsner, 1993 ; Rushton, 1983 ; Ietswaart and Feij, 1989 ) and then by designating these trees as introgressed forms. Our results would rather indicate that the within species variation encompassed these “ intermediate ” forms, despite the interfertility between the two species. Natural hybridization has been shown by mating system analysis with gene markers in natural mixed populations (Bacilieri et al., 1996 ; Kremer et al., 1998). Depending on the marker system used (isozymes or microsatellites) the species and the pollination season, the estimated rate of hybridization varies from 3% to 32 %. Furthermore, interspecific hybridization was also demonstrated by controlled crossings; Kleinschmit and Kleinschmit (2000) reported from 5 to 13 % success in interspecific crossing experiments as compared to 17% to 20% for within species success. However there was no experimental evidence that introgressed forms  exhibited intermediate morphology. In a review of plant hybridization, Rieseberg and Ellstrand (1993) challenged the usual assumption of hybrid intermediacy for morphological characters. In their survey, they noticed that hybrids were a mosaic of phenotypes with parental and intermediate characters rather than just intermediate ones ; they stated that “ from a systematic perspective, the unpredictability of hybrid character expression diminishes the utility of morphological characters for hybrid identification. ”. Observations exploring  these ideas were made in F1 families of interspecific crosses between Q. petraea and Q. robur (Kleinschmit et al., 1995). These authors showed that F1 juvenile hybrids (2 to 5 years old) exhibited leaf morphologies that were similar to the female parent rather than intermediate, whether the female parent was Q. petraea or Q. robur. At least at the juvenile stage, these maternal effects suggested that F1 hybrids did not exhibit intermediate phenotypes.  We lack observations on older material that could have sustained the same conclusions. Maternal effects could decrease as trees get older, and hybrids could exhibit intermediate phenotypes at the mature stage. An alternative interpretation to our results, in comparison to the extensive hybridization that has been observed, was disruptive selection. If there were numerous hybrid seedlings produced, and if local site conditions favored parental phenotypes, hybrids would have been progressively eliminated as the stand grew older. There was no experimental evidence that hybrids were selected against. However there  are numerous reports of a higher proportion of trees with intermediate morphologies at the edges of the natural distribution, particularly under northern latitudes (Cousens, 1963 Olsson, 1975a), and our results confirmed this. If one accepted that intermediate morphologies corresponded to introgressed trees despite the review by Rieseberg and Ellstrand (1993), and that hybrids had higher fitness in new environmental conditions, than Cousens’ and Olsson’s observations could indicate that disruptive selection could vary according to site conditions.
To sum up, our conclusions about the separation of mixed stands in two populations (the two bell shaped curves on Figure 1) corresponding to the two species based on adult leaf morphology, were not in contradiction with the occurrence of natural hybridization. Either hybrids did not persist until the adult stage due to disruptive selection, or, hybrids exhibited parental morphologies even at the mature stage. However, the genetic composition of the two populations that we identified with the multivariate statistical techniques, would have been  quite different according to the inheritance of the leaf morphological traits and the existence or not, of disruptive selection. If there was disruptive selection, whatever the inheritance of leaf morphology traits was, the two populations would have been composed of pure species with no introgressed forms. If morphological traits exhibited predominantly maternal effects and there was no disruptive selection, then the two populations would have been composed of different genetic entities, pure species and introgressed forms. If there was biparental inheritance and no disruptive selection, then the intermediate forms would have been hybrids and there would have been three populations (the two parental ones plus the introgressed forms, the latter being infrequent or indistinguishable from the two others). Further investigations based on a common assessment of molecular markers of species differences and leaf morphology, as well as the study of hybrid fitness in mixed stands, are needed in order to choose among these alternatives.
Acknowledgements  
The study has been carried out with financial support from the Commission of the European Communities, Agriculture and Fisheries (FAIR) specific RTD programme, CT-FAIR 1 PL95-0297, “Synthetic maps of gene diversity and provenance performance for utilization and conservation of oak genetic resources in Europe”. It does not necessarily reflects its views and in no way anticipates the Commission’s future policy in this area.
We are grateful to Cathleen Baldwin, Jan Bovenschen, Fabienne Bourquin, Cedric Demeurie,  Viggo Jensen, Adrain Jordan, Gert Kranenborg, Thomas Küno, Jean Marc Louvet, to the technical research Unit of the Forestry Commisssion and INRA   for assistance in labelling,  mapping trees, and  measuring leaf morphology, to Tom Connolly for statistical advice during the comparative analysis of data sets, and to the  owners of the forests  for giving access to the study stands. 
References

Aas G., 1993 Taxonomical impact of morphological variation in Quercus  robur and Q. petraea : a contribution to the hybrid controversy. Ann. Sci. For. 50 : 107s-114s
Aas G., 1995. Die Behaarung der Blätter von Traubeneiche und Stieleiche (Quercus  petraea und Quercus  robur) : Variabilität und taxonomische Bedeutung. Mitt. Fortsl. Versuchanstalt Rheinland-Pfalz 34 : 297-309
Aas, G., Friedrich K., 1991. Untersuchungen zur morphologischen Unterscheidung von Stiel-und Traubeneicheln. Forstw. Cbl. 110 : 349-357
Bacilieri R., Ducousso A., Kremer A., 1995. Genetical, morphological, ecological and phenological differentiation between Quercus  petraea (Matt.) Liebl. and Q. robur L in a mixed stand of northwest of France. Silvae genetica 44 :1-10
Bacilieri R., Ducousso A., Petit R.J., Kremer A., 1996. Mating system and asymmetric hynridization in a mixed stand of European oaks . Evolution 50 : 900-908
Benzécri J.P., 1992. Correspondence analysis handbook, Dekker M., New York, 665 pp.
Bodénès C., Joandet S., Laigret F., Kremer A., 1997. Detection of genomic regions differentiating two closely related oak species Quercus  petraea (Matt.) Liebl. and Q. robur L.. Heredity 78 : 433-444
Borovics A., 1999. Separability of pedunculate oak and sessile oak aggregate : a contribution to the hybrid and small species estimation. Proceedings of the Forest Research Institute 89: 93-110. 
Carlisle A., Brown A.M.F., 1965. The assessment of the taxonomic status of mixed oaks (Quercus  ssp.) populations. Watsonia 6: 120-127
Cousens J.E. 1963. Variation of some diagnostic characters of the sessile and pedunculate oaks and their hybrids in Scotland. Watsonia 5 : 273-286
Dupouey J.L., and Le Bouler H., 1989. Discrimination morphologique des glands de chênes sessile (Quercus  petraea (Matt.) Liebl) et pédonculé (Q. robur L.). Ann. Sci. For. 46: 187-194
Dupouey J.L., and Badeau V.,  1993. Morphological varaibility of oaks (Quercus  robur L., Quercus  petraea (Matt.) Liebl. and Q. pubescens Willd) in northeastern France : preliminary results. Ann. Sci. For. 50 :35s-40s
Elsner G., 1993. Preliminary variability of oak stands (Quercus  petraea and Quercus  robur) in northern Germany.  Ann. Sci. For. 50 : 228s-232s
Grandjean G., Sigaud P., 1987. Contribution à la taxonomie et à l’écologie des chênes du Berry. Ann. Sci. For. 44 : 35-66
Iestwaart J.H., Feij, A.E. 1989. A multivariate analysis of introgression  between Q. robur and Q. petraea in The Netherlands. Acta Bot. Neerl. 38 :313-325
Kissling P., 1977 Les poils des quatre espèces de chênes du Jura (Quercus  pubescens, Q. petraea, Q. robur et Q. cerris). Ber. Schweiz. Bot. Ges 87 : 1-18
Kissling P., 1980a. Un réseau de corrélations entre les chênes (Quercus L.) du Jura. Ber. Schweiz. Bot. Ges 90 : 1-28 
Kissling P., 1980b. Clef de détermination des chênes médioeuropéens (Quercus  L.) Ber. Schweiz. Bot. Ges 90 : 29-44
Kleinschmit J., and Kleinschmit J.G.R., 2000. Quercus  robur-Quercus  petraea : a critical review of the species concept. Glas, sum. Pokuse 37 : 441-452
Kleinschmit J.R.G., Bacilieri R., Kremer A., Roloff A. 1995. Comparison of morphological traits of pedunculate oak (Q. robur L.) and sessile oak (Q. petraea (Matt.) Liebl. Silvae Genetica 44 : 256-269
Kremer A., Petit R.J., Ducousso A., 1998 Structure of gene diversity, gene flow and gene conservation in Quercus  petraea. In « First Euforgen meeting on social broadleaves », edited by J.Turok, A. Kremer and S. de Vries, IPGRI p 133-144
Minihan V.B.,  Rushton B.S., 1984. The taxonomic status of oaks (Quercus  ssp.) in Breen Wood, Co Antrim, Northern Ireland. Watsonia 15 : 27-32
Moreau F., Kleinschmit J. R.G., Kremer A.. 1994. Molecular differentiation between Q. petraea and Quercus  robur assessed by random amplified DNA fragments. Forest. Genet. 1 :51-64
Muir G., Fleming C.C., Schlötterer Ch., 2000. Species status of hybridizing oaks. Nature 405 : 1016
Olsson U., 1975a. A morphological analysis of phenotypes in populations of Quercus  (Fagaceae) in Sweden. Bot . Not. 128 : 55-68
Olsson U., 1975b. The structure of stellate trichomes and their taxonomic implication in some Quercus  species (Fagaceae). Bot. Not. 128 : 412-424.
 Rieseberg L.H., Ellstrand N.C., 1993. What can molecular and morphological markers tell us about plant hybridization. Critical reviews in Plant Sciences 12 : 213-241
 Rushton B.S., 1978. Quercus  robur L. and Quercus  petraea (Matt.) Liebl. : a multivariate approach to the hybrid problem. 1. Data acquisition, analysis and interpretation. Watsonia : 12 :81-101
 Rushton B.S., 1979. Quercus  robur L. and Quercus  petraea (Matt.) Liebl. : a multivariate approach to the hybrid problem. 2. The geographical distribution of population types. Watsonia : 12 : 209-224
 Rushton B.S., 1983. An analysis of variation of leaf characters in Quercus  robur L. and Quercus  petraea (Matt.) Liebl. population samples from Northern Ireland. Irish Forestry 40 : 52-77
Samuel R., Pinsker W., Ehrendorfer F. 1995. Electrophoretic analysis of genetic variation within and between populations of Quercus  cerris, Q. pubescens, Q. petraea, and Q. robur (Fagaceae) from eastern Austria. Bot Acta 108 : 290-299
Safou O., Saint Martin M., 1989. Le trichome de quelques Quercus  périméditerranéens. Bull. Soc. Bot. Fr., 136 : 291-304
SAS Institute Inc., 1989. SAS/STAT User's Guide, Version 6, SAS Institute Inc., Cary, North Carolina, 4th edit., 2 vol., 1848 p.
Wigston D.L., 1975. The distribution of Q. robur L., Q. petraea (Matt.) Liebl. and their hybrids in southwestern England. 1. The assessment of the taxonomic status of populations from leaf characters. Watsonia 10 : 345-369
Zanetto A., Roussel G., Kremer A., 1994. Geographic variation of inter-specific differentiation between Quercus  robur L. and Quercus  petraea (Matt.) Liebl.. Forest genet. 1 : 111-123

Captions of Figures

Figure 1a

Title
Description of the dimensional leaf morphological traits in Quercus petraea and Q. robur.
Legend
LL: Lamina length  
PL: Petiole length 
LW: Lobe width 
SW: Sinus width 
WP: Length of lamina from base to  widest point
PU : Abaxial laminar pubescence 

Figure 1b

Title
Scoring of the basal shape of the lamina (BS) of Quercus petraea and Q. robur. 

Figure 2

Title
Distribution of the first two synthetic variables according to the different multivariate analyses of leaf morphological traits in Quercus petraea and Q. robur.

Legend
A, C and E correspond to the diagram of the trees along the  first two synthetic variables  (1st variable as x axis, second variable as y axis)
B, D and F represents the distribution of the first synthetic variable 
A and B : Canonical Discriminant Analysis (CDA1 and CDA2 are the  first two synthetic variables of CDA)
C and D : Principal Component Analysis (PCA1 and PCA2 are the  first two synthetic variables of PCA)
E and F : Multiple Correspondence Analysis (MCA1 and MCA2 are the  first two synthetic variables of MCA)

Figure 3

Title
Distribution of the discriminant function values according to the species
Legend
Each tree was assigned to either Q. petraea or Q. robur when the three multivariate methods provided congruent results.  When there was a discrepancy among the methods the tree was considered as unclassified (see text).  

Figure 4

Title
Distribution of the 3025 trees according to petiole length (PL, horizontal axis) and number of intercalary veins (NV, vertical axis). 

Legend
  Q. petraea  (□)
  Q. robur ()
Unclassified (■)

The line represents the zero value of the discriminant function separating the two species.