Local adaptation through genetic differentiation in highly fragmented Tilia cordata populations

Abstract We assessed the level of geographic differentiation of Tilia cordata in Denmark based on tests of 91 trees selected from 12 isolated populations. We used quantitative analysis of spring phenology and population genetic analysis based on SSR markers to infer the likely historical genetic processes within and among populations. High genetic variation within and among populations was observed in spring phenology, which correlated with spring temperatures at the origin of the tested T. cordata trees. The population genetic analysis revealed significant differentiation among the populations, but with no clear sign of isolation by distance. We infer the findings as indications of ongoing fine scale selection in favor of local growth conditions made possible by limited gene flow among the small and fragmented populations. This hypothesis fits well with reports of limited fruiting in the investigated Danish T. cordata populations, while the species is known for its ability to propagate vegetatively by root suckers. Our results suggest that both divergent selection and genetic drift may have played important roles in forming the genetic patterns of T. cordata at its northern distribution limit. However, we also speculate that epigenetic mechanism arising from the original population environment could have created similar patterns in regulating the spring phenology.

. Selection and random drift occur simultaneously in natural populations and have in the case of tree species often resulted in the occurrence of genetically differentiated populations over large geographic ranges (Alberto et al., 2013;Nadeau, Meirmans, Aitken, Ritland, & Isabel, 2016). The balance between neutral processes and natural selection is important, and the ability of a species to adapt to particular set of local growth conditions can be limited if neutral processes dominate the evolution of among population variation (Savolainen, Pyhäjärvi, & Knürr, 2007).

| Genetic differentiation in fragmented landscapes
Habitat fragmentation can have a negative influence on a species' capacity to adapt to new conditions by reducing genetic diversity and increasing inbreeding in small populations (Lowe, Cavers, Boshier, Breed, & Hollingsworth, 2016). However, tree species in fragmented landscapes often maintain connectivity through extensive gene flow in the form of pollen movement and seed dispersal (Breed, Ottewell, Gardner, & Lowe, 2011). On the other hand, a reduced gene flow among populations can be beneficial for local adaptation, because site-specific fitness of adaptive alleles may only lead to local adaptation in the absence of homogenizing gene flow (Sork, 2016) . Trees are long-lived organisms, and it is therefore difficult to estimate effects of divergent selection directly by comparing development over several generations. Instead, phenotyping in common garden trials or provenance field trials are applied to compare the performance of trees that originate from different climatic conditions. In such studies, evidence of natural selection is inferred from the level and pattern of genetic differentiation among origins from divergent climatic conditions (Morgenstern, 1996). More recently, landscape genomics involving correlation of allele frequencies and environmental conditions has been applied to infer effects of selection, although the efficiency of the approach with the present level of genomic data is still discussed (Ćalić, Bussotti, Martínez-García, & Neale, 2015).
Estimation of gene flow across landscapes based on genetic markers is, however, well established, and results from a large number of studies of trees are available (Lowe, Cavers, Boshier, Breed, & Hollingsworth, 2015;Savolainen et al., 2007). A combination of the different approaches is therefore ideal for studying the genetic background of local adaptation of trees in fragmented landscapes (Lepais & Bacles, 2014;Sork et al., 2013).

Tilia cordata
Mill. (small-leaved lime) is native to Europe with its Northern distribution limit at the Southern part of Finland, Sweden, and Norway. Being close to its Northern limits, native Danish T. cordata mostly occur in small isolated populations in ancient forests (Lawesson, 2004). The species was dominating in Danish forests until 2,500 years ago (Hannon, Bradshaw, & Emborg, 2000), and the present fragmented distribution of T. cordata can therefore be a result of relatively early colonization by the species after postglaciation followed by recurring events of local extinction during the recent part of Holocene, a development that may have taken place in other parts of the natural distribution as well (Fineschi, Salvini, Taurchini, Carnevale, & Vendramin, 2003). Gene flow among the remnant populations is likely to be limited at the northern limit (Myking, 2002), especially due to limited production of fertile seeds (Pigott & Huntley, 1981), which is also observed in many of the Danish populations included in this study (Lawesson, 2004). While the pronounced ability of T. cordata to reproduce through root and stump suckers (Koop, 1987a,b) can maintain the population size and genetic diversity within populations, it will result in regeneration without gene flow even among closely located populations.

| Objective of the study
In this study, we investigate the hypothesis that divergent selection (which possibly has occurred during thousands of years) combined with very low level of gene flow among the presently very fragmented populations has led to adaptation at a very local scale in Denmark. Also, we study the level of genetic diversity based on the expectation that the life-history characteristics of the species such as pollination mechanism determine the genetic differentiation happening within and among populations.
In order to separate the effects of genetic drift and selection, we compare the level of population differentiation in the quantitative genetic trait spring phenology (Q ST value) with population differentiation based on putative neutral SSR markers (F ST value).
The presence of divergent natural selection will be indicated by Q ST > F ST and uniform or stabilizing selection by Q ST < F ST , while Q ST = F ST suggests that genetic drift is the major driver behind differentiation among populations in the studied traits (Leinonen, Scott McCairns, O'Hara, & Merilä, 2013;Whitlock & Guillaume, 2009). This is because the Q ST value denotes the entire population differentiation in a given trait due to combined effect of natural selection and neutral processes, while F ST at neutral loci only measures genetic distance between populations arising as a result of neutral processes such as gene flow and drift (De Kort, Vandepitte, & Honnay, 2012). We test if any observed population differentiation in spring phenology can be explained by climatic differences at the locations, or if the patterns mainly reflect geographic distances between populations. We study spring phenology, as there are indications of bud burst advancement in T. cordata (Kramer, 1995), and hence to test if the investigated populations are prone to late spring frost damages in the future. Finally we discuss, what can be learned from the findings in relation to the species' ability to adapt its phenology to the expected changes in the future due to climate change.

| Plant material
The study was based on 91 trees (clones) selected for the Danish genetic resource conservation program (Graudal, Kjaer, & Canger, 1995)

| Quantitative genetic analysis of spring phenology
Bud burst of all trees in the clonal test was assessed on 19 April 2004 and 5 May 2006 using a scale from 1 to 6 based on the development stage of bud burst from fully closed winter buds (score 1) to fully unfolded leaves (score 6). From these bud burst data, population values (least square mean values [LSMeans]) and variance components between populations and between clones within populations were estimated using model (1) given below: where Y ijkl is the bud burst score measured for tree l, μ is the overall mean of the bud burst score, B i is the fixed block effect, P j is the fixed population effect, λ ij is the fixed population by block interaction, C k(j) is the random effect of clones within population, and ε ijkl is the residual. The broad sense heritability (H 2 ) was calculated according to Falconer and Mackay (1996) where V E is the estimated residual (environmental) variance in the clonal trial. We estimated the expected response to selection (R) for bud burst using breeder's equation following Falconer and Mackay (1996) as follows; R = iH 2 √V P, where i is the selection intensity, H 2 is the broad sense heritability, and V P is the phenotypic variance.
We used a scenario based on selection of the 5% and 10% most extreme phenotypes. We compared this measure of expected response from one round of strong selection with the present differences among populations in order to illustrate the magnitude of the present levels of population differentiation in timing of budburst.
The actual level of genetic differentiation in bud burst was calculated following Spitze (1993) as is the variance between populations and V G is the estimated genetic variance of the clones. Q ST values are downwards biased as we use the total genetic (clonal) variance as proxy for the additive genetic variance (Goudet & Büchi, 2006;López-Fanjul, Fernández, & Toro, 2003). The software program ASReml v3.0 was used to estimate (1) Y ijkl = μ+B i + P j +λ ij + C k(j) +ε ijkl , F I G U R E 1 Minimum May temperature and location of the 12 populations (red stars) included in the study (Map provided by Mikael Scharling, Danish Meteorological Institute) and with an outline of four eco-geographic regions. Phenology was assessed in a clonal seed orchard with replications (black star) variance components and provenance values as well as standard errors of broad sense heritability estimates and Q ST estimates (Gilmour et al. 2009).

| Molecular analysis based on SSR markers
Nuclear DNA was extracted from the sampled leaves with the QIAGEN ® DNeasy 96 Plant Kit (Germany), using approximately 40 mg of leaf material and following the manufactures instructions.
The samples were genotyped for nine nuclear microsatellites developed by Phuekvilai and Wolff (2013). Tc23 was not described in Phuekvilai and Wolff (2013), but included in Hansen, Thomsen, and Rasmussen (2014). The microsatellites were amplified in three different multiplex primer mixes using Qiagen multiplex PCR kit (Germany) and following the PCR protocol optimized by Phuekvilai and Wolff (2013). The PCR products were visualized by capillary electrophoresis on an ABI 3130xl sequencer (Applied Biosystems), and size standard GZ500LIZ was applied for reference of fragment size and scored using GeneMapper v.4.0 (Applied Biosystems).
Genetic diversity in the four eco-geographic regions used in the study was estimated for each locus using the following parameters: were calculated in GenAlEx ver. 6.5 (Peakall & Smouse, 2006 while allelic richness was calculated in the HP-Rare 1.1 software (Kalinowski, 2005). Pairwise F ST values between clones and between trees pooled in eco-geographic regions were calculated in GenAlEx 6.502 (Peakall & Smouse, 2005).

| Genetic and phenotypic difference among clones as a function of geographic distance
We used a mantel test implemented in the package ade4 in R version 3.2.2 (Dray & Dufour, 2007)

| Quantitative genetic variation in spring phenology
The maximum distance among the 12 populations studied was only ~300 kilometers, and the difference in their mean annual tempera- p-value = .004) at the 12 original population sites belonging to the four climatically distinct eco-geographic regions in Denmark (Figure 2).
Pairwise geographic distance between populations and their corresponding differences in LSMeans of bud burst score showed a tendency to be related (r = .22; p-value = .06) (Figure 3).

| Molecular genetic variation
Principal coordinate analysis (PCoA) run on GenAlEx revealed no patterns among the 12 populations studied ( Figure S1). Hence, we based our further genetic analysis on the four eco-geographic regions to which the populations belong. A summary of the results from the genetic analysis is given in Table S1. The pairwise F ST values among the clones pooled in the four eco-geographic regions were small, but in general significant (Table 3). Largest F st value (0.037) was found between the southwest and north eco-geographic regions. However, population F st values and their corresponding geographic distance were not correlated (r = .02; p-value = .28) (Figure 4).

| D ISCUSS I ON
The 12 populations in the study were significantly differentiated with respect to spring phenology. Contrary to phenotypic data, can also have been important in shaping the present patterns in this trait (Hutchison & Templeton, 1999;Pickup, Field, Rowell, & Young, 2012). Conclusions are valid since our estimate of Q ST for bud burst is downwards biased, but still much larger than F ST . Common gardens such as the clonal field trials in the present study helps in accounting for possible phenotypic plasticity exhibited by individual genotypes as they are grown in a single environment (Franks et al., 2014;Merilä & Hendry, 2014). However, environmental effects at the original population sites can still induce plastic/epigenetic effects (De Kort et al., 2014;Dewan et al., 2018;Groot, Wagemaker, Ouborg, Verhoeven, & Vergeer, 2018). Such epigenetic differentiation can be population specific and can be involved in adaptation to local environments (Bossdorf, Richards, & Pigliucci, 2008;Herrera & Bazaga, 2010). Hence, it is important to discern between the genetic and epigenetic regulation of adaptive traits among populations before drawing valid conclusions regarding mechanisms behind local adaptation (Grativol, Hemerly, & Ferreira, 2012;Richards, Bossdorf, & Verhoeven, 2010).
In our study, trees from warmer population origins flushed earlier at the test site (regression of population bud burst scores on minimum temperatures in spring at the population sites were significant). The populations grouped within each eco-geographic region showed a clear pattern in their correlation between phenology and temperature in spring. This suggests that these populations are isolated by environment rather than distance (Sexton, Hangartner, & Hoffmann, 2013). In population genetics, it is often observed that the isolation by environment is more important for local adaptation than isolation by distance (Tiffin & Ross-Ibarra, 2014); as under the influence of isolation by environment, it is likely that adaptive alleles undergo divergent selection within populations (Nadeau et al., 2016

| CON CLUS IONS
Our findings exemplifies that population differentiation in trees often occurs as a result of natural selection and neutral processes simultaneously. In our case, the populations have maintained a high level of genetic variation and thereby possess ability to respond to selection based on high level of heritability within populations together with high population differentiation. This suggests the presence of adaptive potential for spring phenology if exposed to strong selection. Plasticity also allows species to adapt to changes in growing conditions rapidly, and we cannot exclude that epigenetics have played an important role in creating the local adaptation (Verhoeven, vonHoldt, & Sork, 2016). Repeated assessments of budburst over years and ideally different climates will be required to compare the importance of genetic variation with phenotypic plasticity. The role of epigenetics is still poorly studied in woody species, but studies in Norway spruce (Picea abies) suggest that in at least some species the mechanism may be very efficient (Yakovlev et al., 2014). In the case of Tiliaceae, no studies have to our knowledge investigated the potential role of epigenetics in local adaptation, and this aspect therefore calls for more research in order to understand the adaptation mechanism in the species and support qualified discussions on the need for interventions as assisted migration or selection in the face of rapid climate change.

ACK N OWLED G M ENTS
The authors would like to thank the Villum Foundation for the financial support to carry out this research as part of the Trees for future forests project (VKR-023063).

AUTH O R CO NTR I B UTI O N
Albin Lobo was responsible for the design of study, analysis and interpretation of data, and writing of the manuscript. Ole Kim Hansen and Eva Ortvald Erichsen were responsible for the marker analysis, interpretation of molecular data, and writing of manuscript. Jon Kehlet Hansen was responsible for supervision of the work, writing and approving of manuscript. Birgitte Jacobsen was responsible for field data collection and statistical analysis of data. Erik Dahl Kjaer was responsible for overall supervision of the study, conception and design of work, data interpretation, writing and approval of the manuscript.

DATA A RCH I V I N G
Data for this study is available at University of Copenhagen-Electronic Research Data Archive (UCPH ERDA).