High Genetic Differentiation among European White Oak Species ( Quercus spp . ) at a Dehydrin Gene

Dehydryn genes are involved in plant response to environmental stress and may be useful to examine functional diversity in relation to adaptive variation. Recently, a dehydrin gene (DHN3) was isolated in Quercus petraea and showed little differentiation between populations of the same species in an altitudinal transect. In the present study, interand intraspecific differentiation patterns in closely related and interfertile oaks were investigated for the first time at the DHN3 locus. A four-oak-species stand (Quercus frainetto Ten., Q. petraea (Matt.) Liebl., Q. pubescens Willd., Q. robur L.) and two populations for each of five white oak species (Q. frainetto Ten., Q. petraea (Matt.) Liebl., Q. pubescens Willd., Q. robur L. and Q. pedunculiflora K. Koch) were analyzed. Three alleles shared by all five oak species were observed. However, only two alleles were present in each population, but with different frequencies according to the species. At population level, all interspecific pairs of populations showed significant differentiation, except for pure Q. robur and Q. pedunculiflora populations. In contrast, no significant differentiation (p > 0.05) was found among conspecific populations. The DHN3 locus proved to be very useful to differentiate Q. frainetto and Q. pubescens from Q. pedunculiflora (FST = 0.914 and 0.660, respectively) and Q. robur (FST = 0.858 and 0.633, respectively). As expected, the lowest level of differentiation was detected between the most closely related species, Q. robur and Q. pedunculiflora (FST = 0.020). Our results suggest that DHN3 can be an important genetic marker for differentiating among European white oak species.


Introduction
The genus Quercus L. comprises roughly 500 species (Aldrich and Cavender-Bares, 2011) subdivided into several distinct monophyletic groups (Manos et al., 1999;Manos and Stanford 2001).Oaks are trees or shrubs which are found in evergreen or deciduous forests throughout the Northern Hemisphere, ranging from temperate to tropical or semi-arid regions (Nixon, 1993).Their high ecological and economic importance is well known (McShea et al., 2007;Aldrich and Cavender-Bares, 2011), being among the most notable forest trees.At the same time, oaks are among the most difficult to distinguish at the species level, mainly due to their phenotypic plasticity and capacity for interspecific hybridization (Gailing and Curtu, 2014).These Quercus features led to debate the biological species concept (Mayr, 1942) and to serve as a model in the development of species ecological concept (Van Valen, 1976).
In the last decade, due to great advances in genetic and molecular analysis, the main issue of differentiation between species has moved at the genome level.Thus, it was proposed that genes can be considered units of species differentiation (Wu, 2001).In this model, the sympatric and interfertile species, as many oaks are, possess a genome with a mosaic of impermeable and permeable regions to gene flow.While the permeable ones share introgressed genes that decrease interspecific differentiation, the impermeable regions are responsible for the maintenance of species integrity because they accumulate divergence in response to selection.This theory is supported by several studies in interfertile plant and animal species (Barton and Gale, 1993;Via and West, 2008), as well as in the Quercus spp.complex (Scotti- Saintagne et al., 2004).From this point of view, one can better understand the antagonistic effects of interspecific gene flow and divergent selection which causes the low genetic differentiation between hybridizing oak species at most genomic regions and the high interspecific differentiation at several 'islands' (outlier regions), respectively.This pattern is consistent with the expected spatial clustering of outlier markers due to ecological speciation (Via, 2009).
was previously outlined at Bejan Forest (Curtu et al., 2007;Curtu et al., 2015) and individuals with intermediate morphology were categorized as the species they most closely resemble.

Genetic analysis
After sample collections in the field, plant material (twigs with buds) was transferred to the laboratory and frozen at -60 ° C. Subsequently, DNA was extracted from buds using the Qiagen DNeasy96 and Plant Mini Kit following manufacturer's instructions and the protocol modified by Toader et al. (2009).The DHN3 marker was amplified using an unlabelled primer pair (5'-TCCATCACTCCCTTCTTCTGA -3' (forward) and 5' -TGTCGCATTACCAAAACCAG -3' (reverse)) and a polymerase chain reaction carried out with a Corbett thermocycler.The PCR program included the following steps: (1) denaturation at 95 °C for 3 min; (2) 30 cycles with a denaturation step at 95 °C for 40 s, primer annealing at 55 °C for 1 min and an elongation step at 72 °C for 2 min; (3) final elongation at 72 °C for 15 min.Subsequently, the PCR products and a 50 bp ladder (Promega) were run in 1.5% agarose gel stained with GelRed (Biotium), by means of horizontal electrophoresis.PCR fragments were run on gel for 10 min at 50 V and 100 min at 90 V. Allele scoring was performed under UV illumination according to the known size of the ladder.

Statistical data analysis
The following summary statistics were estimated within populations and species, at DHN3 locus, using GenAlEx 6.5 (Peakall and Smouse, 2012): number of alleles (Na), number of effective alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He) and Wright's inbreeding coefficient (FIS).In addition, allelic richness (Ar), a diversity measure that is independent of the sample size (El Mousadik and Petit, 1996;Petit et al., 1998) was computed using the program FSTAT 2.9.3.2 (Goudet, 2001), with a rarefaction of n = 41 individuals (the minimum population size found in RO-Pet).After a grouping of populations to corresponding species, differences in levels of Ar, Ho, He and FIS between each pair of species were analyzed using the 'comparison among groups of samples', with 1000 permutations for significance test.Potential deviations from the Hardy-Weinberg equilibrium (HWE) at DHN3 gene were tested for each population using Fisher's global exact tests based on Markov chain analysis (Guo and Thompson, 1992) included in GENEPOP 4.2 (Rousset, 2008).Genotypes variation in each population was calculated in MS Excel.A cluster analysis using the unweighted pair group method with arithmetic mean (UPGMA, Sneath and Sokal, 1973) algorithm was performed based on Nei's (1972) minimum genetic distance.Genetic distances and trees were computed using the POPULATIONS 1.2.32 software (Langella, 1999), with 1000 replications for the bootstrap resampling.Subsequently, the dendrogram of the oak populations was created in MEGA 6 (Tamura et al. 2013).For genetic differentiation among oak populations and species, pairwise FST's for all pairs of populations and species, respectively, were computed using the software ARLEQUIN 3.5.2.2 (Excoffier et al., 2005).The significance of the FST statistics was tested by 10000 permutations.The graphical representations of all pairwise FST were done using a Rfunction (pairFstMatrix.r,Schneider et al., 2000) implemented in ARLEQUIN via Rcmd.With the same software, a The outlier loci, closely linked or even directly under selection (Lind-Riehl et al., 2014), display significantly higher FST values than those expected under selectively neutral conditions, mainly due to different proportion of the alleles among the species (Muir and Schlötterer, 2005;Lind and Gailing, 2013).In addition, the genomic regions containing outliers between Q. petraea and Q. robur indicated lower recombination rates as compared to a control region (Goicoechea et al., 2012).Because such loci might be responsible for adaptive differences among oak species, their identification became very important for proper species assignment.
Recently, the isolation and characterisation of a full-length dehydrin gene (DHN3) and its promoter region in Q. petraea (Vornam et al., 2011) were carried out.The dehydrin genes are putatively involved in plant response to environmental stress and can be useful to examine functional diversity in relation to adaptive variation.The aim of this study was to provide information on the distribution of genetic diversity at DHN3 gene within and among five Quercus species in Romania and to evaluate the potential of this marker for species differentiation in the white oak complex.

Sampling strategy
Two populations for each of the five white oak species (Quercus frainetto Ten., Q. petraea (Matt.)Liebl., Q. pubescens Willd., Q. robur L. and Q. pedunculiflora K. Koch) were selected within their natural distribution in Romania (Table 1).The geographical distance between each pair of species ranged from 120 km (Q.pedunculiflora populations) to 385 km (Q.pubescens populations).In addition, the Bejan Forest, a mixed stand with four white oak species (Q.frainetto, Q. petraea, Q. pubescens and Q. robur) was sampled.In order to avoid, as much as possible, the sampling of related trees, a minimum distance of 50 m between individuals was kept in pure species populations.On the other hand, in the Bejan Forest all trees were sampled on an area of 8,6 ha.Morphological identification was performed following Șofletea and Curtu (2007).Evidence of natural hybridization hierarchical Analysis of Molecular Variance (AMOVA, Excoffier et al., 1992) was employed to examine the partitioning of genetic variation into components: among groups of populations (species in our case), among populations within groups and within populations.Distributions generated from 10000 random permutations were used to estimate the significances of the various variance components.

Genetic variation within species
In 967 oaks genotyped at the DHN3 gene, three alleles shared by all five oak species were found.However, only two alleles were present in each population, but with different frequencies according to the species (Fig. 1).
The fragment's length of the alleles were visually estimated after its relative position to the known size of the ladder lines.Thus, the most frequent alleles were labelled as 274 bp and 310 bp and were found in similar proportions (48.7% and 51.0%, respectively) in the overall population (n = 967 trees).On the other hand, the third allele, labelled as 346 bp, is very rare (0.03% in all samples) and was found only in heterozygous individuals from five populations (Table 2), each one made up of different species.It should be noted that in the initial description of the DHN3 gene in five Q. petraea natural populations, Vornam et al. (2011) found only two alleles: EMBL AM711635 and EMBL AM711636, with PCR fragments length of 656 bp and 620 bp, respectively.Thus, the allele had a 36-nucleotide deletion within the coding region and corresponds to the allele denoted 274 bp in this study.The other allele (656 bp) is indicated as 310 bp in our study.
The allelic richness (Ar), a genetic parameter without population size bias, showed a small variation among populations (Table 2) and no significant difference among species (p > 0.05).On the other hand, higher observed heterozygosity (Ho) and gene diversity (He) values were identified in Q. petraea and Q. pubescens populations as compared to Q. frainetto, Q. robur and Q. pedunculiflora populations (Table 2).
As a consequence, when averaged across populations within each species, Ho differed significantly (p < 0.05) between Q. petraea and Q. frainetto, Q. robur and Q. pedunculiflora, as well as between Q. pubescens and Q. pedunculiflora.The same situation was observed for He, except for the last pair of species.No population showed significant deviations (p < 0.05) from genotypic frequencies expected under the Hardy-Weinberg equilibrium (Table 2), although high values of FIS indicated heterozygote deficit in several populations.Differentiation among species The DHN3 gene showed a high capability to differentiate among four oak species.However, no significant differentiation (p > 0.05) was found among populations of the same species (Fig. 3a).In contrast, slight differentiation in the coding region (FST = 0.015) was reported at the same dehydryn gene among five Q. petraea populations from different altitudes in the French Pyrenees (Vornam et al., 2011).Second, no significant differentiation was found between pure Q. robur and Q. pedunculiflora populations.This result is consistent with the taxonomic relation between these oak taxa, Q. pedunculiflora being a minor species or a subspecies of Q. robur (Curtu et al., 2011b).However, a small (FST = 0.028 and 0.040, respectively) but significant difference (p < 0.05) was identified between the Q. robur population from the mixed stand (BE-Rob) and the Q. pedunculiflora populations (BC-Ped and PU-Ped).Also, all other pairs of populations with different species showed significant genetic differentiation (p < 0.05).In contrast, at species level, all pairwise differentiations were significant (p < 0.05).Q. frainetto and Q. pubescens are markedly differentiated from Q. pedunculiflora (FST = 0.914 and 0.660, respectively) and Q. robur (FST = 0.858 and 0.633, respectively) at DHN3 locus (Fig. 3b).As expected, the lowest level of differentiation was detected between the most closely related species, Q. robur and Q. pedunculiflora (FST = 0.020).However, the latter significant differentiation between this pair of species could be biased by the presence of populations from the mixed stand (Bejan Forest), BE-Rob especially.Indeed, when a complementary differentiation analysis was set up without any population from the Bejan Forest, Q. robur and Q. pedunculiflora could not be significantly differentiated.In addition, the FST values between the remaining pairs of species were similar with those obtained previously, with a slight increase in each case (data not shown).
The high capability of interspecific differentiation at DHN3 gene between almost all pairs of analyzed species suggests that this marker can be an important outlier in the European white oak complex.While only three of the ten pairs of species showed interspecific FST < 0.20 (Fig. 4b), the majority of pairs of species had an interspecific FST > 0.30.In the Bejan Forest, differentiation among four species was high (FST > 0.20) at two nuclear SSRs (ZAG36 and ZAG96) between some pairs of species (Curtu et al., 2007).However, DHN3 performed better than the afforementioned SSRs, but with a similar FST value between Q. petraea and Q. robur (0.33 vs. 0.31 at DHN3 and ZAG96, respectively).Up to date it is known that a low proportion of the oak genome is highly differentiated (Scotti- Saintagne et al., 2004) and few markers can match as outliers.Thus, in a complex genome study with many marker types (AFLPs, SSRs, SNPs, SCARs and isozymes), Scotti- Saintagne et al. (2004) reported only 12% of total loci (47 out of 389) as outliers, with interspecific FST > 0.018 between Q. robur and Q. petraea.Neophytou et al. (2010) found three nuclear SSRs with interspecific FST ranging from 0.18 to 0.38 between the same pair of species.The latter FST value corresponds with the ZAG96 locus, linked with a quantitative trait locus (QTL) associated with petiole length ratio (Saintagne et al., 2004), which showed high interspecific FST values in other oak studies (Muir and Schlötterer, 2005;Curtu et al., 2007).However, a very small FST value (0.03) was reported at the same locus between Q. robur and Q. petraea populations closed to putative refugial regions (Yücedağ and Gailing, 2013).Recently, Guichoux et al. (2013) identified 60 outliers within a set of 262 SNPs which presented interspecific FST values ranging from 0.27 to 0.93 between Q. robur and Q. petraea.These values suggest that the DHN3 gene can also be an outlier among all four species.A similar situation was identified at EST-SSR FIR013 (Lind-Riehl et al., 2014) which displayed a very high interspecific differentiation between Q. rubra and Q. ellipsoidalis (FST = 0.668).This locus is putatively located in COL-1 (CONSTANS-like) gene that could be involved in adaptive divergence and reproductive isolation between these two red oak species.Interestingly, similar allelic variation pattern in DHN3 and COL-1 genes suggest that few alleles have distinct frequencies according to the species.Thus, both genes might be present in genome regions under divergent selection and could be involved in ecological speciation (Nosil, 2012).
In plants, dehydrins belong to group 2 of LEA (late embryogenesis abundant) proteins and are one of the well studied putative dehydration protective molecules.These types of proteins are produced in response to dehydration during  1997), low temperatures or freezing (Welling et al. 2004;Yakovlev et al. 2008), and increased salinity (Close, 1997;Rorat 2006).However, the precise function of these genes remains currently unclear.The DHN3 gene can be classified as a K3 dehydrin at Q. petraea (Vornam et al., 2011) that primarily responds to low temperature (Campbell and Close, 1997).In addition, the Kn-type dehydrins are induced by low temperatures and correlate with plants freezing tolerance capacity (Rorat 2006).Thus, homozygote carriers of short allele (i.e.274 bp in this study) might be less freezing-tolerant than heterozygous or homozyous individuals at longer allele (i.e.310 bp) (Vornam et al., 2011).Expanding the correlation at species level it seems that Q, frainetto and Q. pubescens are less tolerant to freezing than Q. petraea and Q. robur.Indeed, this observation is consistent with the generally accepted ecological requirements of these species.Thus, Q. frainetto and Q. pubescens are found in Southern and Central Europe, while Q. petraea and Q. robur are the most widespread European species, with the latter adapted even to the more continental climate of Northeastern Europe.However, further studies on this gene are required to confirm the assumption that DHN3 gene is involved in freezing tolerance in the Quercus complex.In addition, it should be investigated if dehydrins was involved in the regulation of soluble carbohydrates concentrations, because oak species with higher carbohydrate concentration showed better cold hardiness (Morin et al., 2007).As a consequence of the lack of significant genetic differentiation between Q. robur and Q. pedunculiflora, the grouping of populations for hierarchical AMOVA was based on the results of pairwise FST and were: (1) Q. frainetto, (2) Q. petraea, (3) Q. pubescens, (4) Q. robur and Q.Thus, AMOVA showed that the percentage of molecular variance is almost equally distributed among species (52.6% of total variance) and within populations (47.3% of total variance), while the variance among populations within species is very small and not significant (p = 0.30).Similar results were reported in a gene coding for a non-specific NAD-dependent dehydrogenase (Gömöry, 2000), where the interspecific component of variation accounted for 54.3% of the total gene diversity and 43.5% of total variation was distributed within populations.In two pairs of Q. frainetto and Q. pubescens populations the largest molecular variance is distributed within populations (92.7%) and only 6.2% of the total variance is between species (Curtu et al., 2011a).On the other hand, in a study on three oak species (Q.frainetto, Q. petraea and Q. pubescens), the molecular variance is mostly distributed within the species while a smaller portion of genetic variability was explained by variance among species (86% and 14% of total variance, respectively) (Antonecchia et al., 2015).
The dendrogram constructed using Nei's genetic distances between pairs of populations revealed four main clusters which mostly comprise populations of a single species (Fig. 4), with very high reliability of nodes (100), based on bootstrap resampling.However.Q. pedunculiflora populations were grouped together with Q. robur populations.A relatively small genetic distance was found between Q. frainetto and Q. pubescens.Also, a high genetic affinity was noticed between Q. petraea and Q. robur -Q.pedunculiflora group.Similar results were reported by Curtu et al. (2007) in a mixed stand with four oak species, but based on isozymes and nuclear SSRs as genetic markers.Thus, these convergent results suggest that the grouping of oak species is maintained irrespective of the particular marker type used.

Conclusions
The analysis of the DHN3 gene indicates high levels of genetic differentiation among the four major white oak species occurring in Romania.At the same time Q. robur could not be distinguished from its closest relative, the drought-tolerant Q. pedunculiflora, which is considered a minor or incipient species.Populations cluster together according to species rather than geographic origin except for Q. robur -Q.pedunculiflora complex.The hierarchical analysis of molecular variance and pairwise FST point out that DHN3 is a potential outlier that stands out from the low differentiated genome of European white oak species.Further studies needed to confirm that the DHN3 gene is responsible with freezing tolerance in Quercus complex.

Fig. 1 .
Fig. 1.Allele frequency at DHN3 gene in 10 pure populations and one mixed forest (Bejan) consisting of four oak species.The populations name is displayed at the top of pie chart.The species name is abbreviated as follows: FRA-Q.frainetto, PED-Q.pedunculiflora, PET-Q.petraea, PUB-Q.pubescens, ROB-Q.robur Fig. 2. Allele variation at the DHN3 gene among five oak species in Romania

Table 1 .
Geographic location of the sampled Quercus populations and of the mixed forest of four oak species

Table 2 .
Diversity measures of DHN3 gene within oak populations and species