Geographical Structuring of Quercus robur (L.) Chloroplast DNA Haplotypes in Lithuania: Recolonization, Adaptation, or Overexploitation Effects?

: We studied the maternally inherited chloroplast DNA polymorphism at three microsatellite loci of 157 Quercus robur trees from 38 native populations in Lithuania. We found high diversity of eight haplotypes from the Balkan lineage A (frequency 0.75) and the “German” subbranch of the Balkan lineage A (freq. 0.12), western and eastern Italian lineages C (freq. 0.05 and 0.06, respectively), and Iberian lineage B (freq. 0.03). The haplotypes were geographically well structured (among population differentiation index PhiPT = 0.30, the p -value < 0.001) that is unexpected for such a small territory as Lithuania. We raised a hypothesis on historical overexploitation of oaks by eliminating certain haplotypes in Lithuania, following a drastic felling of oak forests over the last few centuries.


Introduction
Quercus robur (L.) is an important forest key-habitat species, an indicator of the autochthonous forests with immense values for economy, ecosystem services and culture in Europe [1]. In Lithuania alone occupying merely ca. 63 K km 2 , there are 74 settlements bearing the name of "oak" [2]. Revealing the evolutionary history of such an ecologically important species would markedly benefit the basic knowledge on genetic structure of forest tree populations and help to efficiently design genetic biodiversity conservation programs [3]. The efficiency here is understood as capturing rare genetic variants rather than the genetic variation common elsewhere. There is a possibility for these rare variants to reside in divergent evolutionary lineages of forest tree populations [4]. Another issue is the decline of health of Q. robur forests observed over the recent decades in Europe [5][6][7][8][9]. One possible cause of the oak decline could be depleted genetic diversity due to the centuries of exploitation of oak forests in Europe [5].
A strategy to test the oak decline hypothesis is to study richness and structuring of evolutionary lineages of oaks within a certain region [10]. Presence of structuring of evolutionary lineages, that may be difficult to explain due to natural causes, may be a strong indicator of undesirable human interventions such as overexploitation. The evolutionary origin of plants is commonly determined by studying the genetic variation in maternally inherited chloroplast genomes [11]. In oaks, chloroplast DNA is maternally inherited and is passed over generations less affected by stochastic recombination events as in the nuclear genome [11,12]. Mutations in tree plastid genomes are not frequent and the consensus is that certain maternally inherited DNA lineages (haplotypes) are associated to the glacial refugium populations and the postglacial migration lineages of trees expanding since the LGM at 20,000 years BP in late Pleistocene geological epoch [13,14]. Studies on the clear genome [11,12]. Mutations in tree plastid genomes are not frequent and the consensus is that certain maternally inherited DNA lineages (haplotypes) are associated to the glacial refugium populations and the postglacial migration lineages of trees expanding since the LGM at 20,000 years BP in late Pleistocene geological epoch [13,14]. Studies on the fossil pollen abundance in soil sedimental horizons and DNA polymorphism of maternally inherited plastid genomes suggest that genetically distinct zones for postglacial expansion of trees were located in the Iberian and Italian peninsulas, and the Balkans [15,16]. Factors other than the geographic proximity may also affect migration of trees, e.g., river flow direction [17,18]. Based on maternally inherited chloroplast DNA markers the following three main postglacial migration lineages of Q. robur were detected in Europe [14,[19][20][21] (Figure 1): Lineage "A" originating from glacial refugium in the Balkans during the migration accumulated mutations leading to several haplotypes, one of which outbranched this lineage in present-day Germany and is often referred to as "German" haplotype; Lineage "B" out of the Iberian peninsula spreading mainly to France and Britain but also reaching as far as the Scandinavian peninsula; Lineage "C" from southern Italy that mutated into western and eastern subbranches, correspondingly spreading towards north-west and north-east of Europe, the latter being common in the western regions of European Russia [21]. Geophysical barriers and corridors as well as proximity to the glacial refugia were the main determinants of the present-day geographical arrangement of the post-glacial migration lineages of Q. robur. For instance, the corridor between the Alps and the Carpathian Mountain chains as well as the northward alignment of the Carpathians favored the Balkan refugia of plants for migration northwards over the Italian refugia.  [19] (noted by the capital letters A, B, C) and haplotypes as in [14] (noted by capital H) and our study (noted by Haplo). Haplo-4-5 indicates the "German" branch of the Balkan lineage A. The location of Lithuania on the map is highlighted.
The Balkan lineage (A) is the geographically closest migration route into the eastern Baltic region [14] and was found to be more frequent in earlier studies on cpDNA markers  [19] (noted by the capital letters A, B, C) and haplotypes as in [14] (noted by capital H) and our study (noted by Haplo). Haplo-4-5 indicates the "German" branch of the Balkan lineage A. The location of Lithuania on the map is highlighted.
The Balkan lineage (A) is the geographically closest migration route into the eastern Baltic region [14] and was found to be more frequent in earlier studies on cpDNA markers of oaks in the Baltic region [10,22]. These studies reported high chloroplast DNA haplotype diversity with evidence of geographically non-random distribution raising a hypothesis of local adaptation of certain haplogroups. However, ref. [10] focused on broader eastern Baltic gradient with few populations from Lithuania. Both these studies used the PCR_RFLP DNA marker system which is not precise enough to detect certain haplotypes that can be detected with the chloroplast microsatellite markers [20,21,23]. The radiocarbon dating of the fossil pollen in Europe suggests a fast postglacial tree expansion northward at an average rate of 150 km per century (conifers), 100 km/century (Corylus, Ulmus, Alnus) and 40 km/century (Tilia, Quercus, Fagus Carpinus) mainly during 14,000 to 12,000 years BP following the warming at the start of the Holocene geological epoch [15,[24][25][26]. The species composition of European forests at certain periods of time depended on (a) the geographical distance from these glacial migration sources (e.g., [14,27]) and (b) once populated, on the long-term temperature fluctuations usually referred to as climate chronozones affecting the species biological adaptability to certain climates [24]. Therefore, each part of Europe has its own time scale for forest species abundance and composition. Based on range wide radiocarbon dated fossil pollen records [28], oaks populated present-day southern France and Germany during the relatively cooler Preboreal and Boreal climatic chronozones (10,000-8000 years BP) and reached the eastern Baltic region between 8000 and 6000 years BP during the warm Atlantic climatic chronozone. According to [26], the most abundant oak fossil pollen records north of the latitude 47 • N were for 11,000-8000 years BP. However, the oak fossil pollen records in southern Lithuania indicate frequency of oaks up to 1% of the other tree species already during the tundra vegetation period 14,000-13,000 BC [29]. The oldest radiocarbon dated oak fossil wood from mature tree trunks found in the eastern Baltic are of ca. 9000 years BP (fossils conserved at the Baltic seabed or pond peat layer in Birzai, northern Lithuania [30]). Based on fossil pollen, the golden age of oaks in Lithuania was 8000-3000 years BC, with a peak at late Atlantic zone in 6500 years BC [29]. Around 6500 years BC, oaks admixed with other broadleaf trees may have formed several continuous forest tracts and amounted up to 10% of other forest tree species [29]. These large oak dominated forest tracts were located (a) in south-western Lithuania stretching from the central part of Kaunas-Kedainiai via Kalzu Ruda towards Marijampole and southwards, (b) in the western Lithuanian highland (Zemaitija, centered about Rietavas), and (c) in an area from Kaunas stretching eastwards towards Trakai [29]. Since then, following the Subatlantic cooling with expansion of Norway spruce, oak abundance was decreasing gradually down to 5% in the late Subboreal zone (3000 BP) and 1-2% at 500 BP and the present-day.
Strong and dynamic genetic differentiation due to divergent adaptive environments at the range widelevel have led to a pronounced ecotypic structuring within Q. robur species sensu lato [31][32][33]. This strong genetic differentiation [34] provided a vast morphological variation for the taxonomists seeking the sensu stricto taxonomic categorization [35,36]. Nature, however, seems to complicate the hard work of the taxonomists even more by allowing a spontaneous hybridization between the taxonomic units not only within but also among the species [37,38] and not only for Quercus sp. but also for other forest trees [39][40][41]. Nevertheless, there is a possibility for a long-term retention of the morphology markers that are neutral to selection such as some of the leaf or acorn traits in Q. robur [33]. In such case, these neutral marker traits to some degree could be reflected by distinct evolutionary lineages initiated from the divergent glacial refugia [31]. For Q. robur, several intraspecific taxonomic structures stand out: Q. robur subsp. robur (common in central and northern Europe); Q. robur subsp. broteroana (Schwarz) Camus, Q. robur subsp. estremadurensis (Schwarz) Camus, common in Spain (both correspond to the Iberian post-glacial migration lineage); Q. robur subsp. Brutia common in southern Italy (the migration lineages from Italy), Q. robur subsp. slavonica (Gayer) Matyas, concentrated in Slovenia (the migration lineages from the Balkan Peninsula, [42]), Q. robur subsp. Pedunculiflora (Asia Minor, Caucasus, Balkans). Based on historical sources, in the 16th century, Q. robur stands may have occupied 10-20% of the forests in Lithuania [2]. The paleontological records indicate a lower margin of this range [29]. In any case, since the 18th century in the Russian Empire, a severe felling of oaks for shipbuilding and construction reduced the Q. robur stands down to 2-3% of its original range by the 1900s [2]. Therefore, it is likely that humans have contributed to the oak decline phenomena by eroding the genetic diversity when selectively removing oaks of high commercial value and by deforesting areas on rich soils suitable for agriculture, as in south-western Lithuania. Presently, Q. robur stands in Lithuania are small and fragmented and occupy 47.3 thousand ha, making up 2.30% of the remaining forests [42]. Commonly, the commercial value of oaks in Lithuania is poor: forked trees, usually losing the apical dominance at the lower half of the stem, dominate (personal observation). There are initiatives at the level of Ministry of Environment to promote and enlarge areas of Q. robur in Lithuania. The question of a depleted genetic diversity of Q. robur populations in Lithuania is very much relevant, when initiating such large-scale oak restitution initiatives.
The objectives of our study were to assess geographical variation patterns in the maternally inherited chloroplast DNA haplotypes among natural Q. robur populations in Lithuania. To supplement the earlier RFLP-based studies on the evolutionary origin of oaks in the Baltics, here we used a dense population sampling, microsatellite genotyping and discussed the main factors that could have influenced the observed structure of oak haplotypes.

Material and Methods
We sampled leaves (where accessible) or wood from 157 young trees at 38 locations (populations) in natural stands of Q. robur representing all of Lithuania ( Figure 2). The young trees were chosen for easy access to leaves. Per location, we selected ca. 5-10 trees spaced more than 30 m apart. By considering the relatively lower variation at the organelle DNA loci, our approach was to capture a more detailed geographical grid at the cost of lower sample size at each grid point. When interpreting the genetic structure, it is easy to visually merge the grid locations into geographically larger units to obtain higher within unit sample size, if needed. The initial aim was to genotype 10 trees per location, however, due to the genotype scoring quality at most of the locations the sample size is lower (Table 1).  genetic analyzer. The alleles were scored by the aid of a binset in GENEMAPPER soft. v4.1 [45]. We assigned the postglacial migration lineages from [19] to our cpSSR haplotypes by synchronizing our SSR allele scores at the two loci ucd4 and udt4 with those from [20,23] by alignment of the most common alleles as shown in Table 3. In [23] the cpSSR haplotype association with the [19] lineages is defined by analyzing the [19] reference samples in a single plate with the studied samples (Table 3). Figure 2. Geographical distribution of Q. robur chloroplast DNA haplotypes in Lithuania. The haplotypes are explained in the legend, where the haplotype id, allele scores (both as in our study), number of trees after the slash, the haplotypes as in [14], and the post-glacial migration lineages as in [19] are given. The dashed arrows indicate hypothetical migration routes of the postglacial migration linages identified as in [19] assuming natural migration being the main driver of the present distribution of the haplotypes. Table 3. Synchronizing the alleles for the two cpSSR loci from our and the two other studies on Q. robur with the aim to assign the cpSSR haplotypes to the postglacial migration lineages by [19]. Q. robur. The table includes all alleles found at the loci in all three studies.

Reference (Origin of the Material)
Alleles at Locus ucd4, bp Size Shift 1 Alleles at Locus udt4, bp Size Shift 1 Figure 2. Geographical distribution of Q. robur chloroplast DNA haplotypes in Lithuania. The haplotypes are explained in the legend, where the haplotype id, allele scores (both as in our study), number of trees after the slash, the haplotypes as in [14], and the post-glacial migration lineages as in [19] are given. The dashed arrows indicate hypothetical migration routes of the postglacial migration linages identified as in [19] assuming natural migration being the main driver of the present distribution of the haplotypes. We extracted the DNA from leaves or sawdust collected by drilling with an electric bore ca. 1-2 cm into the trunk (bore diameter 0.5 cm). For the DNA extraction, we used a modified ATMAB protocol [43]. For the genotyping we used three chloroplast microsatellite (cpSSR) markers ucd_4, udt_4 and udt_1 [44] ( Table 2). The cpSSR primers were selected based on the finding that the pair of markers ucd4 and udt4 is sufficient for distinguishing between the main postglacial migration lineages detectable by PCR-RFLP by [13,20]. The allele number (Na) and the fragment size are given as found in our study. Ta is the primer annealing temperature (TD means touchdown used). All primers from [45].
The DNA concentration taken for the PCR was adjusted to approximately 30-50 ng/µL. The PCR amplification was carried in a 15 µL reaction mix containing 7.5 µL of the Plat-inum™ Multiplex PCR Master Mix (ThermoFisher Scientific, Waltham, MA, USA), 3 µL of RNAse free water, 1.5 µL of Primer Mix (2 µM each primer), 1 µL of DNA, 1 µL of PVP (1%) and 1 µL of BSA 20 mg/mL. The PCR was run with Applied Biosystems (Waltham, MA, USA) thermo cycler GeneAmp PCA System 9700 as follows: an initial pre-denaturation step at 95 • C for 15 min, followed by 26 cycles of denaturation at 95 • C for 30 s, annealing at 57 • C for 1 min 30 s, and extension at 72 • C for 30 s, followed by the final extension step at 60 • C for 30 min. Failed or unclear amplification was repeated. The fragments were separated by performing a capillary electrophoresis on the ABI PRISM™ 310 genetic analyzer. The alleles were scored by the aid of a binset in GENEMAPPER soft. v4.1 [45].
We assigned the postglacial migration lineages from [19] to our cpSSR haplotypes by synchronizing our SSR allele scores at the two loci ucd4 and udt4 with those from [20,23] by alignment of the most common alleles as shown in Table 3. In [23] the cpSSR haplotype association with the [19] lineages is defined by analyzing the [19] reference samples in a single plate with the studied samples (Table 3). Table 3. Synchronizing the alleles for the two cpSSR loci from our and the two other studies on Q. robur with the aim to assign the cpSSR haplotypes to the postglacial migration lineages by [19]. Q. robur. The table includes all alleles found at the loci in all three studies.  1 Size shift is accounted for the alleles in our study as the reference (+1 bp means that an allele is larger for 1 bp in comparison to the allele in our study). * most common allele at each locus in the three studies.

Reference (Origin of the Material) Alleles at Locus ucd4, bp
We constructed the multilocus haplotypes and calculated the haplotype frequency for each population. The genetic differentiation among the populations was tested by AMOVA based on genetic distances with GENALEX v.6.5 [46]. The haplotype geographical Forests 2021, 12, 831 7 of 12 distribution was visualized by displaying the haplotype proportions for each population on the map of Lithuania. To analyze the geographical structure based on the allele frequencies, we used the Monmonier's algorithm allowing for establishing barriers along a significant shift in the allele frequency within a landscape implemented in soft. BARRIER v.2 [47]. The programme (a) creates a Delaunay triangulation plot between the sampled populations, (b) calculates genetic distances (DA genetic distance [48] in our case) associated with each edge in the plot, and (c) creates growing barriers along the largest genetic distances on the plot, the barriers are ranked based on the magnitude of the differentiation.

Results
All three cpSSR loci were polymorphic with two to four alleles yielding eight different haplotypes (Table 4). Based on the alignment for haplotype identity, we found the representatives of all three main postglacial migration lineages A, B and C in Lithuania (Figure 2). Most of the haplotypes (74.5%) originated from the Balkan refugium (coded as Haplo-3 in our study, corresponding to lineage A [19] or haplotypes H10-12 as in [14]). The haplotypes Haplo-3 and Haplo-4 formed the second most abundant haplogroup with 11.5% frequency. This haplogroup originates from the Balkans but outbranched from the main lineage by a mutation in central Europe (likely present-day Germany [20]). The haplotype Haplo-6 with 6.4% frequency represents the eastern Italian lineage common in Belarus and western parts of European Russia (the lineage C, H2). The western Italian haplogroup (the lineage C, H1, in our study the haplotypes Haplo-7 and Haplo-8) occurred at 4.4% frequency. Haplo-1 was found in a single individual in ROKI population but was absent in [20,23] (Table 4). The AMOVA revealed very strong population genetic differentiation based on the cpSSR markers with the PhiPT value as high as 0.3 (the p-value < 0.001) and population variance component of 30% of the total variation in the molecular data (Table 5). The geographical structuring of the chloroplast DNA haplotypes in the country clearly was not random (Figure 2). The most frequent Balkan lineage A (Haplo-3 in our study) was found in most of the populations with relatively higher proportions in the southern part of the country (Haplo-3 in white, Figure 2). The "German" branch of the Balkan lineage A (Haplo-4-5 in our study) spread mainly into northern Lithuania only. The eastern Italian lineage C was found exclusively in the southeastern part of the country (Haplo-6 in black, Figure 2). The western Italian branch of the C lineage, however, was exclusively present in north-eastern Lithuania (Haplo-7 and Haplo-8 in blue diagonals, Figure 2). Finally, in a single seaside population of JUOD, we found four individuals originating from as far as the Iberian lineage B (Haplo-2 in yellow, Figure 2). According to the BARRIER analysis, the following main geographically consistent gradients for significant allele frequency shifts were observed (in the order of decreasing significance): (1) separating the southeastern Lithuanian Q. robur populations (likely due to the presence of the eastern Italian lineage C), (2) separating western Lithuania (likely due to western Italian C and Iberian B lineages) and (3) separating north-eastern Lithuania (likely due to "German" sub-branch of Balkan lineage A) (Figure 2).

Discussion
We used the most frequent allele at each locus to synchronize the oak cpSSR haplotypes with those in [20,23] (Table 3). It is a subjective method, however, the following arguments indicate that we may be near to the correct associations: (a) the most frequent Balkan linage A in our study was also most frequent in [20,23], (b) the number alleles in the ucd4 and udt4 loci was the same as in [23], and (c) for the capillary electrophoresis and the allele scoring, we used the same hardware and software as in [20,23].
As in the previous studies on the Baltic oaks [10,22], we found all three main evolutionary lineages of oaks in Lithuania, among which the Balkan lineage A was the most frequent, and the geographical distribution of the haplotypes was in good agreement with the two above mentioned studies. This lineage diversity confirms repeated oak recolonization events into the eastern Baltic Sea region as discussed by [10]. The national gene conservation program could be adjusted to capture each lineage by at least one genetic reserve.
The geographical barriers of significant shifts in the maternally inherited Q. robur chloroplast DNA allele frequencies among Q. robur populations in Lithuania ( Figure 3) correspond well to the geographical arrangement of the haplotypes ( Figure 2). Basically, these barriers outlined the western Italian lineage C in western Lithuania, the "German" sub-branch of Balkan lineage A in north-east and eastern Italian lineage C in the southeast of Lithuania by leaving central Lithuania for the dominance of the Balkan lineage A (Figure 2).
The geographical distribution of the maternally inherited Q. robur haplotypes in Lithuania clearly has a pattern. This pattern raises questions on why in the light of large geographical scale of evolutionary migration events into such small territory (ca. 65 K km 2 ), there is such discrete geographical structuring of the maternally inherited chloroplast DNA haplotypes? The first question is why the Balkan lineage A is the most widespread and frequent in Lithuania? The Balkan lineage A was also the most frequent cpDNA haplotype in Lithuania based on the PCR-RLFP genotyping by [22]. It is likely that owing to the shortest geographical distance and perhaps favorability of the northward migration flow of major rivers such like Wisla, Order, and Nemunas, the Balkan refugium lineage had an advantage compared to its nearest competitors from the Italian lineage C. In addition, the geographical positioning of the Alps is less favorable for northwards migration from the Italian refugia than that of Carpathians for migration from the Balkan refugia. In agreement with [22], we found the eastern Italian Lineage C in eastern Lithuania. This result may be interpreted as a relatively later north-westwards migration of the eastern Italian lineage from the south-eastern European plains of Russia [21]. It is likely that the eastern Italian lineage C had recently reached the Baltic region, because it is rare elsewhere westwards in Lithuania ( Figure 2). This finding may help to interpret postglacial migration rates.
zation events into the eastern Baltic Sea region as discussed by [10]. The national gene conservation program could be adjusted to capture each lineage by at least one genetic reserve.
The geographical barriers of significant shifts in the maternally inherited Q. robur chloroplast DNA allele frequencies among Q. robur populations in Lithuania (Figure 3) correspond well to the geographical arrangement of the haplotypes (Figure 2). Basically, these barriers outlined the western Italian lineage C in western Lithuania, the "German" sub-branch of Balkan lineage A in north-east and eastern Italian lineage C in the southeast of Lithuania by leaving central Lithuania for the dominance of the Balkan lineage A (Figure 2).  We found the haplotypes of the Iberian lineage B in four trees only, all located in the seaside spit of Neringa. The other studies reported the Iberian haplotypes of Q. robur in north-western Poland (close to the southern connection of the seaside spit to the continent) and eastern Sweden [10,19]. Thus, the samples of the lineage B we found in Neringa could be a continued natural migration from the Polish populations northwards, or less likely bird mediated travelers over the Baltic sea from Sweden. Additionally, there is a possibility of artificial introduction of oaks of the lineage B, because the seaside spit of Neringa formerly was a part of Eastern Prussia, where introduction of forest trees from other parts of the German state were common [49].
The geographical distribution of the remaining haplotypes, however, is less easy to explain by natural processes of migration alone. The first unexpected result under natural migration is a very compact concentration of the western Italian lineage C in the western highlands of Zemaitija. If considering natural migration of this haplotype northwards from Poland via south-western Lithuania, it should have entered Lithuania and spread over a much broader range from southern to western Lithuania. However, with a comparably dense sampling grid we were able to detect this western Italian lineage only in the very north-western part of the country. The second, seemingly odd geographical distribution concerns the "German" sub-branch of the Balkan lineage A. Based on the earlier studies [14,20], this sub-branch of the Balkan lineage arose somewhere in present-day Germany and we assume that it migrated from Germany north-eastwards into Lithuania. If considering a natural migration pathway, as with the western Italian lineage, we should be able to find more representatives of the "German" sub-branch of the Balkan lineage in the southern and central parts of Lithuania. However, this was not the case in our study.
What could be the main reasons for such a geographically biased distribution of the western Italian and the "German" lineages in Lithuania? An explanation could be different forestry practices and seed sources used in the former Prussian and Russian imperial states [50], as is the case of Slavonian oak introduced into North Rhine-Westphalia of Germany in the second half of the 19th century [20]. However, this is unlikely to be the main reason because the former East Prussian province occupied the very seaside and south-western parts of Lithuania, far shorter distances than both the western Italian and the "German" lineages were located in the present study. Another more feasible explanation could be commercial exploitation of oak trees with long, straight, and clear boles that intensified in Lithuania since the 16th century and on to meet the continuing demand of quality oaks for multiple purposes [2]. In support of this hypothesis [51], reported reduced differentiation among Q. robur populations in the region of western Europe under stronger commercial exploitation of forests. Stronger historical forest felling happened in south-western Lithuania [2], where population differentiation is relatively lower and Balkan lineage A dominates (Figure 2). This hypothesis could be true, given the presence of associations between the stem morphotype and the chloroplast DNA haplotype of oak trees. If true, then over centuries, certain commercially fit haplotypes may have gradually been eliminated in southwestern and central Lithuania owing to better trade conditions and stronger domestication of land than in the other more forested regions of the country.
Differential adaptation of the Q. robur haplotypes to eco-climatic gradients of Lithuania is also possible, as discussed in [22]. However, the surface of Lithuania is flat with the altitudinal rage from ca. 10 m.a.s.l. at the seaside to ca. 250 m.a.s.l. on the highlands, the longitudinal and latitudinal gradients are small, spanning ca. 350-400 km. Such gradients apparently are too narrow to cause such a discrete adaptation-based geographical structuring of cpSSR haplotypes as observed in our study. Here again, a study on the haplotype-morphotype associations in Q. robur would be of value.
In conclusion, we found high diversity of Q. robur chloroplast DNA haplotypes for such a small territory as Lithuania containing haplogroups from the four major European lineages: Iberian, western and eastern Italy, and the Balkans. The haplotypes were geographically well structured which is difficult to explain by natural migration alone. This finding is in favor of the hypothesis of commercial elimination of certain haplotypes over centuries of exploitation of oak forests in the region. Repeated migration events may have also contributed to the present-day haplotype structure. For further testing this oak overexploitation hypothesis, a study on haplotype-morphotype associations of oaks would be of value.