# MOLECULAR FOOTPRINTS OF LOCAL ADAPTATION IN TWO MEDITERRANEAN CONIFERS

Delphine Grivet <fnref n="1" />; Federico Sebastiani <fnref n="2" />,<fnref n="3" />; Ricardo Alía <fnref n="1" />,<fnref n="5" />;
Thomas Bataillon <fnref n="4" />; Sara Torre <fnref n="3" />; Mario Zabal-Aguirre <fnref n="1" />;
Giovanni G. Vendramin <fnref n="3" />; Santiago C. González-Martínez <fnref n="1" />,<fnref n="5" />

**Corresponding Author:**

Santiago C. González-Martínez

Department of Ecology & Genetics, Center of Forest Research, CIFOR-INIA, Carretera de la Coruña km 7.5, E28040 Madrid, Spain

e-mail: santiago@inia.es

This study combines neutrality tests and environmental correlations to identify non-neutral patterns of evolution in candidate genes related to drought stress in two closely-related Mediterranean conifers, *Pinus pinaster* Ait. and *Pinus halepensis* Mill. Based on previous studies, we selected 12 amplicons covering six candidate genes that were sequenced in a large sample spanning the full range of these two species. Neutrality tests relatively robust to demography (*DHEW* compound test and *ML-HKA* test) were used to detect selection events at different temporal scales. Environmental associations between variation at candidate genes and climatic variables were also examined. These combined approaches detected distinct genes that may be targeted by selection, most of them specific to only one of the two conifers, despite their recent divergence (< 10 Ma). An exception was *4-coumarate: CoA ligase (4cl)*, a gene involved in the production of various important secondary products that appeared to play a role in local adaptation processes of both pines. Another remarkable result was that all significant environmental correlations involved temperature indices, highlighting the importance of this climatic factor as a selective driver on Mediterranean pines. The ability to detect natural selection at the DNA sequence level depends on the nature and the strength of the selection events, on the timescale at which they occurred and on the sensitivity of the methods to other evolutionary forces that can mimic selection (e.g. demography, population structure). Using complementary approaches can help to capture different aspects of the evolutionary processes that govern molecular variation at both intra- and interspecific levels.

**Key words:** Natural selection, neutrality tests, environmental associations, candidate genes, drought stress, 4cl.

Published in *Molecular Biology and Evolution* (2011) 28: 101-116, with doi:10.1093/molbev/msq190

Tables S1, S2, S3, S4, S5, S6, S7 and S8, and Figures S1, S2 and S3, along with the corresponding supporting information legends and references, can be downloaded here: http://mbe.oxfordjournals.org/content/28/1/101/suppl/DC1

## Introduction

Identifying candidate genes underlying genetic differences for adaptive traits can help to understand how species have adapted to their environment and to predict how they will respond to future climatic changes, which is a special concern in regions such as the Mediterranean Basin, where a substantial decrease in precipitation and a pronounced warming is expected in the near future (Petit, Hampe, and Cheddadi 2005; Giorgi and Lionello 2008). The ability to detect the footprints of natural selection in population DNA sequence samples depends on the nature and the strength of the selection events (Nielsen 2005), on the evolutionary scale at which they occur (Zhai, Nielsen, and Slatkin 2009) and on the sensitivity of the methods to other evolutionary forces that can mimic selection (e.g. demography, population structure; Biswas and Akey 2006).

Positive selection, which drives the increase in frequency of advantageous mutations, is of particular interest because it underlies local adaptation. Recent or on-going positive selection events can be detected using methods based on polymorphism within species (Biswas and Akey 2006; Zhai, Nielsen, and Slatkin 2009). Recently, powerful methods based on site- and haplotype- frequency spectra that are relatively insensitive to demography or population structure have been developed (Zeng, Shi, and Wu 2007). Detecting selection events over a wider evolutionary scale (e.g. recurrent positive selection or ancient selective sweeps), can be efficiently done using methods contrasting within species polymorphism with among species divergence, as these methods are less sensitive to

<footnote n="1">¹ Department of Forest Ecology & Genetics, Center of Forest Research, CIFOR-INIA, Carretera de la Coruña km 7.5, E28040 Madrid, Spain.</footnote>

<footnote n="2">² Department of Agricultural Biotechnology, Genexpress, University of Florence, Via della Lastruccia 14, I50019 Sesto Fiorentino (FI), Italy.</footnote>

<footnote n="3">³ Plant Genetics Institute, Division of Florence, National Research Council, via Madonna del Piano 10, I50019 Sesto Fiorentino (FI), Italy.</footnote>

<footnote n="4">⁴ BIRC - Bioinformatics Research Center, Aarhus University, C.F. Møllers Allé 8, DK8000 Aarhus C, Denmark</footnote>

<footnote n="5">⁵ Sustainable Forest Management Research Institute, University of Valladolid-INIA.</footnote>
assumptions regarding recombination or demography. The *HKA* (Hudson-Kreitman-Aguadé) neutrality test (Hudson, Kreitman, and Aguade 1987) and its extensions (e.g. Wright and Charlesworth 2004) are among the most powerful tests to detect positive selection (Zhai, Nielsen, and Slatkin 2009). In addition, the action of natural selection can be reflected in statistical associations between genetic and environmental data (Manel *et al.* 2003; Hancock *et al.* 2008; Coop *et al.* 2009). Several studies have shown that environmental heterogeneity influences the distribution of genetic diversity across plant populations (see Savolainen, Pyhjäjärvi, and Knürr 2007, for forest trees; Nakazato, Bogonovich, and Moyle 2008; Montesinos *et al.* 2009). Detecting clinal genetic variation along environmental gradients can thus provide evidence for the action of selection (Gram and Sork 2001; Vasemagi and Primmer 2005; Parisod and Christin 2008). Approaches that detect adaptive clinal variation are very attractive because they can directly link environmental gradients with genotypes and phenotypes. However, detecting environmental associations depends on the spatial scale examined and the associations have to be controlled for selectively neutral processes that can also generate clines (e.g. Storz 2002; Hancock *et al.* 2008; Coop *et al.* 2009; Keller *et al.* 2009). Using environmental association approaches in conjunction with neutrality tests helps capturing different aspects of the evolutionary processes that govern molecular variation.

In the present study, we use the combined strategy outlined above to investigate patterns of polymorphism

for a set of candidate genes related to drought stress in maritime pine (also known as cluster pine, *Pinus pinaster* Ait.), and in Aleppo pine (*Pinus halepensis* Mill.). *P. pinaster* is an economically important pine native to the occidental Mediterranean Basin and the European Atlantic front (Figure 1). This species is ecologically versatile, growing in a variety of substrates, in a wide range of elevations and under a range of Mediterranean and Atlantic climate regimes (semi-arid to very humid; Table 1). *P. halepensis*, a pine closely related to maritime pine, is one of the most abundant and widespread pine species of the Mediterranean Basin (Figure 1). This species, which also has a wide ecological breadth, is well adapted to dry conditions as well as high-intensity fire regimes (Tapias *et al.* 2004), which makes it a species of choice for afforestation in xeric Mediterranean areas. These two conifers, because they have extensive and partially overlapping distributions across the Mediterranean Basin, may possibly display adaptation to similar sets of contrasted environmental conditions, i.e. regional climate types and spatial variability. Among the key environmental factors to Pinus distribution in Mediterranean climates, temperature and rainfall have been shown to constitute good biological descriptors (Richardson and Bond 1991; Soto *et al.* 2010). Additionally, several common garden studies have documented differences in drought tolerance and/or gene expression among Mediterranean pine provenances (for *P. pinaster* see Nguyen and Lamant 1989; Costa *et al.* 1998; Chambel, Climent, and Alía 2007; Aranda *et al.* 2009; for *P. halepensis* see Atzmon, Moshe, and Schiller 2004; Sathyan, Newton, and Loopstra 2005;

<figure><img src="image_63.png" /><figcaption>FIGURE 1</figcaption></figure>

Location of the 12 *P. pinaster* and the 9 *P. halepensis* populations and the distribution of the two species. Grey areas represent the two groups (Western

Mediterranean and North African) defined by chloroplast SSRs (see text for details).
Voltas *et al.* 2008). Overall, these studies point to Mediterranean pines as an interesting system to study local adaptation mediated by abiotic stress response. In that context, studying allelic (SNP-allele and haplotype) frequency in drought-related candidate genes along temperature or precipitation clines can help to improve our understanding of pine adaptation to environmental heterogeneity.

Although *P. pinaster* and *P. halepensis* share various ecological requirements and may have responded to their environment in a similar manner, they present quite contrasted biogeographic and demographic histories that can also have a strong effect on the pattern of variation of the candidate genes examined. *Pinus pinaster* has a long-presence in the western Mediterranean Basin surviving the last glaciations in multiple refugia located in southeastern Spain, northern Africa and the Atlantic coast of Portugal (Bucci *et al.* 2007). In contrast, *P. halepensis* would have undergone long-range colonization of the western Mediterranean from ancient populations located in Greece and Turkey (Bucci *et al.* 1998; Morgante, Felice, and Vendramin 1998; Grivet *et al.* 2009). In addition, the

dynamics of *P. halepensis* populations throughout the Mediterranean Basin highly depend on forest fires. This species-specific information on population dynamics has to be taken into account when looking at the molecular footprint of selection.

Based on previous population genomic studies identifying drought tolerance genes potentially under selection (see Eveno *et al.* 2008; Grivet *et al.* 2009) as well as on gene expression studies identifying genes affected by drought stress (Watkinson *et al.* 2003; Rani *et al.* 2009), we examine here patterns of nucleotide variation in an ortholog set of targeted candidate genes in *P. pinaster* and *P. halepensis*. We sample a broad geographic range of populations to: i) provide new evidence about positive selection acting on these genes, ii) identify the time scale at which selection events may have occurred, and iii) examine which environmental factors underlie molecular signatures of selection in both species. Our pluralistic approach provides insights on the adaptive strategy of two conifers that live under the same Mediterranean climate, but present distinct demographic, (re)colonization and life histories.

<caption>TABLE 1</caption>

Spatial coordinates and climatic variables for *P. pinaster* and *P. halepensis* populations

<table><thead><tr><th rowspan="2">Population</th><th rowspan="2">Country</th><th colspan="3">Spatial variables</th><th colspan="8">Climatic variables</th></tr><tr><th>Altitude (m)</th><th>Latitude (degrees)</th><th>Longitude (degrees)</th><th>AMT (°C)</th><th>TS</th><th>MTWM (°C)</th><th>MTCM (°C)</th><th>AP (mm)</th><th>PWM (mm)</th><th>PDM (mm)</th><th>PS</th></tr></thead><tbody><tr><td colspan="13"><strong>P. pinaster</strong></td></tr><tr><td>Arenas de San Pedro</td><td>Spain</td><td>733</td><td>40.194822</td><td>-5.116213</td><td>14.2</td><td>668.9</td><td>33.4</td><td>1.4</td><td>1318</td><td>199</td><td>12</td><td>60.8</td></tr><tr><td>Cazorla</td><td>Spain</td><td>1100</td><td>37.919675</td><td>-2.925765</td><td>11.5</td><td>650.2</td><td>30.6</td><td>-1.4</td><td>1257</td><td>179</td><td>11</td><td>58.3</td></tr><tr><td>Coca</td><td>Spain</td><td>800</td><td>41.254705</td><td>-4.497827</td><td>12.3</td><td>655.9</td><td>31.2</td><td>-0.6</td><td>454</td><td>55</td><td>15</td><td>30.0</td></tr><tr><td>Cómpeta</td><td>Spain</td><td>903</td><td>36.834265</td><td>-3.953989</td><td>14.0</td><td>512.1</td><td>27.5</td><td>4.7</td><td>899</td><td>132</td><td>5</td><td>64.1</td></tr><tr><td>Olba</td><td>Spain</td><td>1002</td><td>40.173309</td><td>-0.622966</td><td>12.4</td><td>600.7</td><td>28.3</td><td>0.6</td><td>509</td><td>63</td><td>24</td><td>31.5</td></tr><tr><td>Oria</td><td>Spain</td><td>1223</td><td>37.531165</td><td>-2.351138</td><td>13.1</td><td>633.1</td><td>30.6</td><td>0.4</td><td>357</td><td>46</td><td>5</td><td>46.2</td></tr><tr><td>Quatretonda</td><td>Spain</td><td>435</td><td>38.971645</td><td>-0.358844</td><td>15.3</td><td>553.7</td><td>30.3</td><td>3.8</td><td>777</td><td>120</td><td>9</td><td>52.1</td></tr><tr><td>Mimizan</td><td>France</td><td>19</td><td>44.134167</td><td>-1.303056</td><td>13.2</td><td>490.9</td><td>24.9</td><td>3.1</td><td>1235</td><td>149</td><td>62</td><td>24.3</td></tr><tr><td>Pinia</td><td>France</td><td>15</td><td>42.021083</td><td>9.464861</td><td>15.6</td><td>515.5</td><td>27.0</td><td>6.3</td><td>585</td><td>80</td><td>8</td><td>49.6</td></tr><tr><td>Tabarka</td><td>Tunisia</td><td>121</td><td>36.958397</td><td>8.703792</td><td>17.7</td><td>573.4</td><td>31.8</td><td>6.7</td><td>916</td><td>162</td><td>4</td><td>72.3</td></tr><tr><td>Tamrabta</td><td>Morocco</td><td>1758</td><td>33.600000</td><td>-5.016667</td><td>11.8</td><td>665.3</td><td>31.3</td><td>-3.2</td><td>721</td><td>96</td><td>8</td><td>54.4</td></tr><tr><td>Sidi Meskour</td><td>Morocco</td><td>1931</td><td>31.439375</td><td>-6.903864</td><td>11.4</td><td>711.3</td><td>32.5</td><td>-5.2</td><td>514</td><td>71</td><td>5</td><td>56.4</td></tr><tr><td colspan="13"><strong>P. halepensis</strong></td></tr><tr><td>Cabanellas</td><td>Spain</td><td>210</td><td>42.235556</td><td>2.790000</td><td>14.8</td><td>575.5</td><td>27.4</td><td>3.9</td><td>713</td><td>97</td><td>35</td><td>26.0</td></tr><tr><td>Carratraca</td><td>Spain</td><td>650</td><td>36.841111</td><td>-4.834444</td><td>15.7</td><td>568.6</td><td>30.6</td><td>4.4</td><td>693</td><td>106</td><td>3</td><td>69.4</td></tr><tr><td>S'avall</td><td>Spain</td><td>10</td><td>39.287222</td><td>3.047778</td><td>16.8</td><td>525.9</td><td>28.7</td><td>6.8</td><td>566</td><td>92</td><td>5</td><td>54.7</td></tr><tr><td>Tarrasa</td><td>Spain</td><td>117</td><td>41.466667</td><td>2.100000</td><td>15.9</td><td>551.1</td><td>27.7</td><td>5.8</td><td>619</td><td>83</td><td>29</td><td>29.5</td></tr><tr><td>Imperia</td><td>Italy</td><td>109</td><td>43.900000</td><td>8.050000</td><td>15.1</td><td>550.5</td><td>27.5</td><td>4.6</td><td>804</td><td>107</td><td>18</td><td>38.5</td></tr><tr><td>Zaouia Ifrane</td><td>Morocco</td><td>1512</td><td>33.570000</td><td>-5.140000</td><td>11.5</td><td>668.7</td><td>31.2</td><td>-3.8</td><td>849</td><td>123</td><td>8</td><td>57.4</td></tr><tr><td>Aures Beni Melloul</td><td>Algeria</td><td>936</td><td>35.166667</td><td>6.833333</td><td>13.4</td><td>657.2</td><td>31.2</td><td>-0.1</td><td>375</td><td>53</td><td>10</td><td>36.1</td></tr><tr><td>Elea</td><td>Greece</td><td>155</td><td>37.766667</td><td>21.533333</td><td>16.7</td><td>592.0</td><td>30.8</td><td>5.4</td><td>808</td><td>150</td><td>6</td><td>77.1</td></tr><tr><td>Shaharia</td><td>Israel</td><td>236</td><td>31.600000</td><td>34.833333</td><td>20.0</td><td>496.4</td><td>31.8</td><td>8.0</td><td>397</td><td>100</td><td>0</td><td>116</td></tr></tbody></table>

AMT: annual mean temperature; TS: temperature seasonality (STD *100); MTWM: maximum temperature of the warmest month; MTCM: minimum temperature of the coldest month;

AP: annual precipitation; PWM: precipitation of the wettest month; PDM: precipitation of the driest month; PS: precipitation seasonality (CV *100).
# Materials and Methods

## Study species

*Pinus pinaster* Ait.: Maritime pine populations have been assessed using various molecular markers: chloroplast and mitochondrial (e.g. Burban and Petit 2003; Bucci *et al.* 2007), as well as nuclear (e.g. Salvador *et al.* 2000; Eveno *et al.* 2008). In particular, chloroplast markers were able to identify different gene clusters related to the history of the species (Bucci *et al.* 2007). Because of its economic importance, various programs of genetic improvement have been developed in *Pinus pinaster*, in particular for the Atlantic provenances. Within this framework, several adaptive traits, such as growth, tolerance to drought and cold, and resistance to pests and diseases, have been the subject of genetic variability studies. Numerous approaches have been used to unravel the basis of quantitative traits in maritime pine: genetic (genetic cartography), physiologic (mechanisms implicated in traits), functional and structural genomic (candidate genes and proteins) and population genetic and genomic (genes under selection).

*Pinus halepensis* Mill.: Aleppo pine genetic variability has been studied with both chloroplast and nuclear markers (Bucci *et al.* 1998; Morgante, Felice, and Vendramin 1998; Grivet *et al.* 2009). A recent study revealed that the pattern of polymorphism observed in some candidate genes related to drought-tolerance in this species reflected long-range colonization and possibly natural selection during range expansion (Grivet *et al.* 2009). Various common garden studies have assessed Aleppo

pine intraspecific variability in order to study the role of ecological factors in shaping adaptive strategies. Some experiments revealed adaptive variation to climate (total precipitation and dry season duration; Voltas *et al.* 2008) or reproductive features (Climent *et al.* 2008), highlighting thus the selective role of climate variables in determining population and family fitness in this species.

## Sampling

Twelve populations of *P. pinaster* (77-122 individuals, depending on the gene; Table S1) and nine populations of *P. halepensis* (72-93 individuals; Table S1) were collected spanning the full range of each species (Figure 1, Table 1). Populations were selected considering not only spatial distribution but also environmental heterogeneity in both species, prioritizing populations that represent contrasted environments (Table 1). Our sampling covers also different soil (siliceous, calcareous) and Mediterranean forest types. Finally, we included representations of the five traditional varieties or landraces described in maritime pine (Resch 1974): *mesogeensis* (e.g. Olba), *atlantica* (e.g. Mimizan), *corteensis* (e.g. Pinia), *maghrebiana* (e.g. Tamrabta) and *renoui* (e.g. Tabarka), as well as all genecological groups normally considered in Aleppo pine: Western Europe, Eastern Europe and Northern Africa.

From each population, cones were collected from mother trees separated by at least 50 m without any phenotypic selection. Seeds from each mother tree were kept in individual paper bags and stored at 4°C in a dry environment till DNA extraction (see below).

<caption>TABLE 2</caption>

Gene diversity (all sites) for 12 amplicons from six putative candidate genes across the twelve *P. pinaster* populations and the nine *P. halepensis* populations

<table><thead><tr><th rowspan="2">Amplicon</th><th colspan="7">P. pinaster</th><th colspan="7">P. halepensis</th></tr><tr><th>N</th><th>L</th><th>S</th><th>Θ<sub>π</sub></th><th>Θ<sub>w</sub></th><th>K</th><th>He</th><th>N</th><th>L</th><th>S</th><th>Θ<sub>π</sub></th><th>Θ<sub>w</sub></th><th>K</th><th>He</th></tr></thead><tbody><tr><td><i>Ip31-Pt</i></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>122</td><td>456</td><td>12</td><td>6.26</td><td>4.89</td><td>6</td><td>0.717</td><td>79</td><td>353</td><td>3</td><td>2.48</td><td>1.72</td><td>4</td><td>0.460</td></tr><tr><td>b</td><td>111</td><td>560</td><td>7</td><td>3.91</td><td>2.37</td><td>3</td><td>0.446</td><td>74</td><td>488</td><td>7</td><td>3.93</td><td>2.94</td><td>4</td><td>0.464</td></tr><tr><td><i>Ip33-Pp</i></td><td>97</td><td>449</td><td>9</td><td>4.29</td><td>3.89</td><td>9</td><td>0.798</td><td>90</td><td>375</td><td>1</td><td>0.17</td><td>0.53</td><td>2</td><td>0.065</td></tr><tr><td><i>dhn2-Pp</i></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>120</td><td>472</td><td>12</td><td>6.06</td><td>4.74</td><td>11</td><td>0.807</td><td>89</td><td>448</td><td>2</td><td>0.20</td><td>0.88</td><td>3</td><td>0.088</td></tr><tr><td>b</td><td>85</td><td>596</td><td>21</td><td>8.42</td><td>7.03</td><td>11</td><td>0.811</td><td>93</td><td>346</td><td>4</td><td>0.55</td><td>2.26</td><td>5</td><td>0.105</td></tr><tr><td><i>dhn2-Ps</i></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>92</td><td>743</td><td>27</td><td>8.10</td><td>7.13</td><td>14</td><td>0.899</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td>b</td><td>80</td><td>513</td><td>15</td><td>5.92</td><td>5.90</td><td>14</td><td>0.855</td><td>92</td><td>457</td><td>11</td><td>1.93</td><td>4.73</td><td>6</td><td>0.204</td></tr><tr><td><i>dhn5-Ps</i></td><td>77</td><td>449</td><td>4</td><td>4.00</td><td>1.81</td><td>5</td><td>0.735</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td><i>4cl-Pt</i></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>91</td><td>551</td><td>13</td><td>5.67</td><td>4.64</td><td>4</td><td>0.678</td><td>88</td><td>461</td><td>10</td><td>5.88</td><td>4.30</td><td>13</td><td>0.805</td></tr><tr><td>b</td><td>92</td><td>256</td><td>7</td><td>3.56</td><td>5.37</td><td>5</td><td>0.291</td><td>79</td><td>688</td><td>21</td><td>10.67</td><td>6.18</td><td>13</td><td>0.793</td></tr><tr><td>c</td><td>87</td><td>543</td><td>4</td><td>0.69</td><td>1.46</td><td>4</td><td>0.341</td><td>78</td><td>528</td><td>3</td><td>2.47</td><td>1.15</td><td>5</td><td>0.697</td></tr><tr><td>d</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>72</td><td>452</td><td>3</td><td>2.80</td><td>1.37</td><td>4</td><td>0.484</td></tr><tr><td><i>Total</i></td><td>1054</td><td>5588</td><td>131</td><td>-</td><td>-</td><td>86</td><td>-</td><td>834</td><td>4596</td><td>65</td><td>-</td><td>-</td><td>59</td><td></td></tr><tr><td><i>Mean</i></td><td>-</td><td>-</td><td>-</td><td>5.17</td><td>4.48</td><td>-</td><td>0.671</td><td>-</td><td>-</td><td>-</td><td>3.11</td><td>2.61</td><td>-</td><td>0.416</td></tr></tbody></table>

N: number of sequences; L: total analyzed length in base pair; S: polymorphic sites; Θ<sub>w</sub>: Tajima's nucleotide diversity × 10⁻³ (Tajima 1989) and Θ<sub>π</sub>: Watterson's nucleotide diversity (Watterson 1975) per site (× 10⁻³); K: number of haplotypes; He: haplotype diversity; na: amplicons that failed to transfer across species.
## Spatial, environmental and genetic data

**Spatial, environmental and genetic data.** Three spatial variables were recorded for each population: altitude, latitude and longitude (Table 1). Eight climatic variables were considered: annual mean temperature (AMT), temperature seasonality (expressed as the standard deviation across months multiplied by 100; TS), maximum temperature of the warmest month (MTWM), minimum temperature of the coldest month (MTCM), annual precipitation (AP), precipitation of wettest month (PWM), precipitation of driest month (PDM), and precipitation seasonality (coefficient of variation; PS). Climatic data for *P. pinaster*'s Iberian populations were obtained from a functional phytoclimatic model based on raw data from meteorological stations (Gonzalo 2007). Climatic data for *P. pinaster*'s non-Iberian populations as well as for all *P. halepensis* populations were obtained from the WorldClim - Global Climate Data at 5 minute resolution (Hijmans et al. 2005) (Table 1). Graphical pairwise correlations between these 11 spatial and environmental variables are presented in Figure S1.

**DNA extraction and candidate gene sequencing.** Genomic DNA from *P. pinaster* (haploid) megagametophytes was extracted with a modified Dellaporta *et al*. (1983) protocol. DNA extraction for *P. halepensis* was carried out as reported in Grivet *et al*. (2009).

Candidate genes related to drought tolerance were originally identified on the basis of functional studies performed in *Pinus taeda* and other conifers, or derived from model species such as *Arabidopsis thaliana*, as described elsewhere (González-Martínez *et al.* 2006; Eveno *et al.* 2008; Grivet *et al.* 2009; Wachowiak, Balk, and Savolainen 2009). Altogether the candidate genes selected for this study belong to three well-known and relatively small multigene families: the ASR (*Ip31-Pt* and *Ip33-Pt*), dehydrin (*dhn2-Pp*, *dhn2-Ps* and *dhn5-Ps*) and 4-coumarate: CoA ligase (*4cl-Pt*) families. The ASR family is named after the Abscisic acid (ABA), stress and ripening response (Frankel *et al.* 2006). These proteins are also induced by water deficit stress (WDS). The dehydrin gene, *dhn2-Pp*, previously described in maritime pine (Eveno *et al.* 2008), is not orthologous to the one described in Scots pine, i.e. *dhn2-Ps* (Wachowiak *et al.* 2009).

Previously published primer pairs of these putative candidate genes were transferred from *P. tadea* (Pt), *P. sylvestris* (Ps) or *P. pinaster* (Pp) to either *P. pinaster* or *P. halepensis* or both, and three new primer pairs, which overlap with pre-described amplicons, were also designed in order to extend the coverage of these target genes (Table S2). PCR conditions for the two pines studied here are given in Table S3. Outgroup sequences for each amplicon were obtained from GenBank: *P. taeda* for amplicons *Ip31-Pt*, *Ip33-Pt*, *dhn2-Pp* and *4cl-Pt*, and *P. sylvestris* for *dhn5-Ps* and one *dhn2-Ps* amplicon; or produced with our newly designed primers (*P. nigra* for

the second *dhn2-Ps* amplicon) (Table S2). Sequence data from this article have been deposited with the EMBL/GenBank Data Libraries under accession nos. HM481479-HM483364.

PCR products were purified and sequenced from both ends following standard protocols (for *P. pinaster*, see Eveno *et al.* 2008; for *P. halepensis*, see Grivet *et al.* 2009). Multiple sequence alignments and manual adjustments were done using SeqMan v7 (DNASTAR Lasergene software) and BIOEDIT.
(http://www.mbio.ncsu.edu/BioEdit/page2.html).

**Chloroplast microsatellites.** To control for associations between candidate genes and environmental data that could be due to neutral processes, we included variation of chloroplast microsatellites (the only neutral genetic markers available for the two species at the wide range scale) as a covariate in the multivariate logistic regressions (see *Environmental associations* section below). We performed population-based Principal Component Analysis (PCA) on chloroplast markers, and kept the

FIGURE 2

Predicted and observed patterns of clinal variation at 4cl-Pt_c. Sigmoid curves denote the clines in either allelic (SNP position given) or haplotypic (haplotypes represented by capital letters) frequencies under a logistic regression model. The logistic regression model uses a climatic variable (MTCM or TS) after controlling for neutral processes (i.e. using the PCs of neutral molecular markers as covariates) for *P. pinaster* and *P. halepensis*. MTCM: minimum temperature of the coldest month; TS: temperature seasonality.

<figure><img src="image_64.png" /></figure>
<caption>TABLE 3a</caption>

Polymorphism, divergence and neutrality tests (silent sites) in *P. pinaster* for the phylogeographical groups defined with chloroplast microsatellites (see text for details).

<table><thead><tr><th>Amplicon</th><th>S</th><th>D</th><th>H<sub>norm</sub></th><th>DHEW P-value</th><th>ML-HKA's K (P. taeda)</th><th>K<sub>s</sub> (P. taeda)</th><th>ML-HKA's K (P. halepensis)</th><th>K<sub>s</sub> (P. halepensis)</th></tr></thead><tbody><tr><td colspan="9"><strong>Western Mediterranean a</strong></td></tr><tr><td>Ip31-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>8</td><td>2.172</td><td>-0.635</td><td>0.958</td><td>0.961</td><td>38.350</td><td>1.292</td><td>20.310</td></tr><tr><td>b</td><td>6</td><td>1.834</td><td>-1.380</td><td>0.913</td><td>0.698</td><td>44.250</td><td>0.720</td><td>24.360</td></tr><tr><td>Ip33-Pp</td><td>6</td><td>1.133</td><td>0.911</td><td>0.984</td><td>0.696</td><td>60.060</td><td>1.839</td><td>15.990</td></tr><tr><td>dhn2-Pp</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>7</td><td>0.269</td><td>-0.853</td><td>0.406</td><td>1.605</td><td>38.160</td><td>1.307</td><td>31.580</td></tr><tr><td>b</td><td>19</td><td>0.186</td><td>0.825</td><td>0.727</td><td>3.256 *</td><td>31.940</td><td>1.975</td><td>39.720</td></tr><tr><td>dhn2-Ps</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>15</td><td>0.631</td><td>0.703</td><td>0.562</td><td>ptna</td><td>ptna</td><td>phna</td><td>phna</td></tr><tr><td>b</td><td>8</td><td>-0.126</td><td>0.232</td><td>0.209</td><td>ptna</td><td>ptna</td><td>d</td><td>d</td></tr><tr><td>dhn5-Ps</td><td>2</td><td>1.960</td><td>0.907</td><td>ptna</td><td>ptna</td><td>phna</td><td>phna</td><td>phna</td></tr><tr><td>4cl-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>2</td><td>0.422</td><td>0.067</td><td>0.610</td><td>0.307</td><td>31.480</td><td>0.250</td><td>26.870</td></tr><tr><td>b</td><td>2</td><td>-1.041</td><td>-3.173 *</td><td>0.014 *</td><td>0.741</td><td>43.930</td><td>0.446</td><td>41.220</td></tr><tr><td>c</td><td>1</td><td>0.214</td><td>0.565</td><td>0.663</td><td>0.542</td><td>23.080</td><td>0.000</td><td>1.190</td></tr><tr><td>d</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td colspan="9"><strong>North Africa b</strong></td></tr><tr><td>Ip31-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>7</td><td>-0.276</td><td>-3.065 *</td><td>0.336</td><td>0.878</td><td>39.300</td><td>0.952</td><td>22.280</td></tr><tr><td>b</td><td>6</td><td>-2.046 *</td><td>-3.813 **</td><td>0.073</td><td>0.820</td><td>44.250</td><td>0.880</td><td>26.150</td></tr><tr><td>Ip33-Pp</td><td>5</td><td>0.014</td><td>0.926</td><td>0.899</td><td>0.614</td><td>57.860</td><td>1.946</td><td>13.040</td></tr><tr><td>dhn2-Pp</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>4</td><td>0.799</td><td>-0.554</td><td>0.615</td><td>0.868</td><td>36.610</td><td>0.792</td><td>28.620</td></tr><tr><td>b</td><td>13</td><td>1.095</td><td>0.361</td><td>0.718</td><td>2.036</td><td>32.780</td><td>1.362</td><td>33.600</td></tr><tr><td>dhn2-Ps</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>6</td><td>0.190</td><td>-0.232</td><td>0.322</td><td>ptna</td><td>ptna</td><td>phna</td><td>phna</td></tr><tr><td>b</td><td>6</td><td>-1.116</td><td>-1.223</td><td>0.085</td><td>ptna</td><td>ptna</td><td>d</td><td>d</td></tr><tr><td>dhn5-Ps</td><td>2</td><td>0.688</td><td>-0.809</td><td>0.540</td><td>ptna</td><td>ptna</td><td>phna</td><td>phna</td></tr><tr><td>4cl-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>8</td><td>2.797</td><td>0.411</td><td>0.995</td><td>1.279</td><td>39.160</td><td>0.442</td><td>28.880</td></tr><tr><td>b</td><td>3</td><td>2.266</td><td>0.090</td><td>0.982</td><td>1.303</td><td>48.190</td><td>1.017</td><td>34.730</td></tr><tr><td>c</td><td>2</td><td>0.201</td><td>0.387</td><td>0.494</td><td>1.039</td><td>25.240</td><td>0.000</td><td>3.410</td></tr><tr><td>d</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr></tbody></table>

D: Tajima's D (Tajima 1989); H<sub>norm</sub>: Fay and Wu's normalized H (Zeng et al. 2006); DHEW test: compound test (Zeng et al. 2007); ML-HKA's K: Selection parameter of the Maximum-Likelihood multilocus Hudson-Kreitman-Aguadé neutrality tests (Wright and Charlesworth 2004); the species used as outgroup is given in parenthesis; S: number of segregating sites; K<sub>s</sub>: average proportion of nucleotide differences between species per silent site ×10<sup>-3</sup>.

na: amplicons that failed to transfer across species; ptna: *P. taeda* outgroup sequence not available; phna: *P. halepensis* outgroup sequence not available; d: discarded because possible paralog.

P-values [*P* neutrality test (neutral) ≤ *P* neutrality test (observed)]: 0 < *** < 0.001 < ** < 0.01 < * < 0.05

a Western Mediterranean: Arenas de San Pedro, Cazoria, Coca, Compeña and Olba.

b North Africa: Tabarka, Tamrabta and Sidi Meskour.

Principal Components (PCs) that explained the majority of the overall inertia of the data. For *P. pinaster*, the three first PCs of the 16 most common haplotypes (determined on the basis of five chloroplast microsatellites) – accounting for 74% of the overall variation – were extracted from the lattice dataset (430 records × 16 variables) of Bucci *et al.* (2007). For populations that did not match any of the 430 locations in this lattice, the mean of the closest 3-4 locations was computed (Table S4). For *P. halepensis*, PCs were computed in R (R Development Core Team) for each of the nine populations based on the 22 haplotypes determined from three chloroplast microsatellites (our unpublished data) and the three first PCs, which explain 82% of the overall variation, were used for subsequent analysis (Table S4).

Chloroplast microsatellites were also used to determine groups with similar evolutionary history to conduct neutrality tests (Figure 1). A previous study based on chloroplast microsatellites pointed at eight gene pools in *P. pinaster* (Bucci *et al.* 2007). Within the sampling used in the present study, Mimizan (Continental French lineage) and Pinia (Northern Italy and Corsican lineages) represent single lineages. Tabarka (Tunisia), Tamrabta (Moroccan Middle Atlas) and Sidi Meskour (Moroccan High Atlas) are all part of a North African lineage that combines western and eastern African origins. The rest of *P. pinaster* populations sampled are part of a wide Western Mediterranean lineage, except for Quatretonda and Oria that are considered marginal populations and may thus present a distinct evolutionary history (González-Martínez *et al.* 2007a; Eveno *et al.* 2008). In *P. halepensis*, eastern
<caption>TABLE 3b</caption>

Polymorphism, divergence and neutrality tests (silent sites) in *P. halepensis* for the phylogeographical groups defined with chloroplast microsatellites (see text for details)

<table><thead><tr><th>Amplicon</th><th>S</th><th>D</th><th>H<sub>norm</sub></th><th>DHEW P-value</th><th>ML-HKA's K (P. taeda)</th><th>K<sub>s</sub> (P. taeda)</th><th>ML-HKA's K (P. halepensis)</th><th>K<sub>s</sub> (P. halepensis)</th></tr></thead><tbody><tr><td colspan="9"><strong>Western Mediterranean a</strong></td></tr><tr><td>lp31-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>2</td><td>1.689</td><td>-1.714</td><td>0.860</td><td>0.597</td><td>34.590</td><td>0.452</td><td>20.960</td></tr><tr><td>b</td><td>4</td><td>0.535</td><td>-0.948</td><td>0.501</td><td>0.889</td><td>48.280</td><td>1.380</td><td>23.090</td></tr><tr><td>lp33-Pp</td><td>1</td><td>-0.629</td><td>-4.343 *</td><td>0.270</td><td>0.145 *</td><td>69.100</td><td>0.361</td><td>14.020</td></tr><tr><td colspan="9"><strong>dhn2-Pp</strong></td></tr><tr><td>a</td><td>0</td><td>nps</td><td>nps</td><td>nps</td><td>nps</td><td>55.770</td><td>nps</td><td>44.020</td></tr><tr><td>b</td><td>0</td><td>nps</td><td>nps</td><td>nps</td><td>nps</td><td>57.180</td><td>nps</td><td>38.220</td></tr><tr><td colspan="9"><strong>dhn2-Ps</strong></td></tr><tr><td>a</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td>b</td><td>2</td><td>-0.595</td><td>0.317</td><td>0.382</td><td>ptna</td><td>ptna</td><td>d</td><td>d</td></tr><tr><td>dhn5-Ps</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td colspan="9"><strong>4cl-Pt</strong></td></tr><tr><td>a</td><td>6</td><td>0.752</td><td>-1.222</td><td>0.599</td><td>2.236</td><td>43.260</td><td>1.618</td><td>22.500</td></tr><tr><td>b</td><td>16</td><td>2.386</td><td>-0.320</td><td>0.985</td><td>2.279</td><td>52.000</td><td>1.242</td><td>30.660</td></tr><tr><td>c</td><td>2</td><td>1.701</td><td>-1.305</td><td>0.867</td><td>1.350</td><td>31.970</td><td>2.976</td><td>9.440</td></tr><tr><td>d</td><td>3</td><td>2.324</td><td>-0.981</td><td>0.965</td><td>0.674</td><td>54.780</td><td>ppna</td><td>ppna</td></tr><tr><td colspan="9"><strong>North Africa b</strong></td></tr><tr><td>lp31-Pt</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a</td><td>2</td><td>1.818</td><td>-1.053</td><td>0.898</td><td>0.819</td><td>34.140</td><td>1.123</td><td>12.710</td></tr><tr><td>b</td><td>2</td><td>1.882</td><td>0.602</td><td>0.909</td><td>0.530</td><td>49.000</td><td>0.924</td><td>22.700</td></tr><tr><td>lp33-Pp</td><td>0</td><td>nps</td><td>nps</td><td>nps</td><td>nps</td><td>69.410</td><td>nps</td><td>9.550</td></tr><tr><td colspan="9"><strong>dhn2-Pp</strong></td></tr><tr><td>a</td><td>1</td><td>-0.774</td><td>0.289</td><td>0.313</td><td>0.329</td><td>56.280</td><td>0.525</td><td>19.310</td></tr><tr><td>b</td><td>1</td><td>-1.147</td><td>0.160</td><td>0.189</td><td>0.215</td><td>57.320</td><td>0.205</td><td>33.600</td></tr><tr><td colspan="9"><strong>dhn2-Ps</strong></td></tr><tr><td>a</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td>b</td><td>0</td><td>nps</td><td>nps</td><td>nps</td><td>ptna</td><td>ptna</td><td>d</td><td>d</td></tr><tr><td>dhn5-Ps</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td colspan="9"><strong>4cl-Pt</strong></td></tr><tr><td>a</td><td>4</td><td>1.306</td><td>-0.201</td><td>0.791</td><td>2.072</td><td>44.030</td><td>1.815</td><td>22.450</td></tr><tr><td>b</td><td>14</td><td>2.316</td><td>0.044</td><td>0.981</td><td>3.251 *</td><td>50.760</td><td>1.951</td><td>31.290</td></tr><tr><td>c</td><td>2</td><td>1.077</td><td>0.221</td><td>0.834</td><td>1.665</td><td>27.720</td><td>3.309</td><td>5.200</td></tr><tr><td>d</td><td>3</td><td>-1.494</td><td>-6.584 ***</td><td>0.009 **</td><td>0.860</td><td>60.020</td><td>ppna</td><td>ppna</td></tr></tbody></table>

D: Tajima's D (Tajima 1989); H<sub>norm</sub>: Fay and Wu's normalized H (Zeng et al. 2006); DHEW test: compound test (Zeng et al. 2007); ML-HKA's K: Selection parameter of the Maximum-Likelihood multilocus Hudson-Kreitman-Aguadé neutrality tests (Wright and Charlesworth 2004); the species used as outgroup is given in parenthesis; S: number of segregating sites; K<sub>s</sub>: average proportion of nucleotide differences between species per silent site × 10<sup>-3</sup>.

na: amplicons that failed to transfer across species; ptna: *P. taeda* outgroup sequence not available; ppna: *P. pinaster* outgroup sequence not available; d: discarded because possible paralog; nps: monomorphic fragment.

P-values (*P* neutrality test (neutral) ≤ *P* neutrality test (observed)): 0 < *** < 0.001 < ** < 0.01 < * < 0.05

Mediterranean populations (Elea in Greece and Shaharia in Israel) present a higher level of cpSSR diversity and substantial genetic differentiation from the rest and are considered single-population lineages (Grivet *et al.* 2009); Zaouia Ifrane (Morocco) and Aures Beni Melloul (Algeria) form part of the North African group, while all four Spanish populations and the Italian population of Imperia belong to the Western Mediterranean group. These regions (North Africa and Western Mediterranean) define relatively homogeneous zones of the Mediterranean Basin in terms of soil and climate (Barbéro *et al.* 1998).

## Statistical analyses

*Gene diversity and divergence.* Number of segregating sites (S), nucleotide diversity statistics (θπ, Tajima 1989; θw, Watterson 1975), number of haplotypes (K), and

haplotypic diversity (H<sub>e</sub>) were computed for both North African and Western Mediterranean groups using scripts kindly provided by S.E. Ramos-Onsins (Department of Genetics, Faculty of Biology, University of Barcelona, Spain) and the program DnaSP v5 (Librado and Rozas 2009). Average divergence per site (K<sub>all</sub>) for each geographical group was computed between each of the pines and the outgroup *Pinus taeda*, as well as between *P. pinaster* and *P. halepensis* themselves, using scripts also provided by S.E. Ramos-Onsins.

*Neutrality tests.* All neutrality tests were performed at the regional level (see group definition above) considering all sequenced gametes in each amplicon, using one sequence from another pine species as outgroup when appropriate, and considering all substitutions as well as only silent sites (except for the *MK* test where all sites must be used, see

<footnote>ª Western Mediterranean: Cabanellas, Carratraca, Imperia, S'avall and Tarrasa.</footnote>

<footnote>ᵇ North Africa: Aures Beni Melloul and Zaouia Ifrane.</footnote>
below). First, we conducted neutrality tests based on within-species population genetic data with (to determine ancestral states) or without outgroup, including Tajima's D test (Tajima 1989), Fay and Wu's *H* normalized test (Zeng *et al.* 2006), both based on the site-frequency spectrum (SFS), Ewens-Watterson *EW* homozygosity test (Watterson 1978), and the *DHEW* neutrality test (Zeng, Shi, and Wu 2007) - a compound test that combines the properties of *D*, *H* normalized and *EW* tests - using scripts provided by K. Zeng (State Key Laboratory of Biocontrol and Key Laboratory of Gene Engineering of the Ministry of Education, Sun Yat-sen University, Guangzhou, China). Second, we used methods that compare patterns of polymorphism within species and divergence between

species. The McDonald-Kreitman (MK) test, which compares the ratio of nonsynonymous and synonymous mutations between and within species (McDonald and Kreitman 1991), was computed using DnaSP v5 (Librado and Rozas 2009). ML-HKA (Maximum-Likelihood multilocus Hudson-Kreitman-Aguadé) neutrality tests were computed using the mlhka program (Wright and Charlesworth 2004). This test is an extension of the HKA test that allows to test specific locus or groups of loci against a group of neutrally evolving genes for their levels of polymorphism and divergence. First, each amplicon was compared to all the rest and the selection parameter *K* was computed. The *K* parameter measures the degree to which diversity is increased or decreased with respect to

FIGURE 3

**Upper section:** *dhn5-Ps* haplotype distribution showing the proportion of haplotypes *C* (in grey), *E* (in black) and all others (in white) in *P. pinaster*. The size of the pie is proportional to the sample size at each location.

**Lower section:** expected (sigmoid curves) and observed covariation in haplotype frequencies of *dhn5-Ps_C* and *dhn5-Ps_E* and climatic variables.

Expected curves were obtained under a logistic regression model including a single environmental variable (ATM or MTWM) after controlling for neutral processes (i.e. using the PCs of neutral molecular markers as covariates) for *P. pinaster*. AMT: annual mean temperature; MTWM: maximum temperature of the warmest month.

<figure><img src="image_65.png" /></figure>
divergence. Second, for significant tests, the analysis was
repeated considering only amplicons from different genes
in order to avoid potential linkage. All runs consisted in
100,000 cycles of the Markov chain; long runs of
1,000,000 cycles were also computed for significant tests
as well as different starting values for the divergence time
parameter (T). However, no differences were found in the
results and only basic runs are presented here. Standard
likelihood ratio tests (LRTs) were used to compare among
models considering one amplicon under selection at a
time and the null model of no selection (degrees of
freedom equal one). The ML-HKA test was performed
using different outgroups to gain insights on evolutionary
times at which positive selective wave(s) may have taken
place: i) with *P. taeda* as outgroup for both *P. pinaster* and
*P. halepensis* amplicons. In this case, the time frame would
correspond to the split between the *Trifoliae* section (to
which *P. taeda* belongs) and the *Pinus* section (to which
the Mediterranean pines belong) in the Oligocene (-25 Ma;
Gernandt *et al.* 2005; Gernandt *et al.* 2008). Three
amplicons were discarded from these analyses because
they had *P. sylvestris* or *P. nigra* as outgroups (two from
dhn2-Ps, and one from dhn5-Ps); ii) with *P. halepensis* as
outgroup for *P. pinaster* amplicons and *P. pinaster* as
outgroup for *P. halepensis* amplicons. The time frame, in
this case, would correspond to the diversification of
Mediterranean pines (in the Miocene, -10 Ma; Gernandt *et
al.* 2008). One amplicon was discarded from this series of
tests (dhn2-Ps_b; see Table S2 for amplicon nomenclature)
due to bad alignment between species and, thus, the
possible amplification of a paralogous fragment in *P.
halepensis* (dhn2-Ps_b from *P. halepensis* had more
similarity with the outgroup, *P. nigra*, than with its close
relative *P. pinaster*).

*Environmental associations.* We tested if SNP-allele or haplotype frequencies at candidate loci correlated with climatic variables using logistic regression. We first carried out series of univariate logistic regressions to test for association between SNP/haplotype frequencies and environmental variables using the program SAM (Joost, Kalbermatten, and Bonin 2008). We considered a correlation as significant only when two likelihood ratio tests (G and Wald tests) rejected the null hypothesis of no association between the genetic and the environmental variables (at the 5% level). A strict Bonferroni correction was applied to correct for multiple testing of univariate models.

In a second step, whenever a SNP/haplotype was found
significantly correlated with an environmental variable, we
performed a multivariate logistic regression using neutral
marker PCs (see above) as covariates (together with the
SNP and the environmental variables) to control for neutral
processes (e.g. postglacial migrations, geographical
isolation) that may have also generated clines in the
absence of any local adaptation. Multiple logistic
regressions were performed in *R* (R Development Core

Team) using the *glm* function (assuming that allelic counts
at SNPs were binomially distributed) for each
SNP/haplotype separately. Model selection was carried out
using the *drop1* function and likelihood ratio tests (LRTs).
Whenever the models compared were not nested (and
LRTs were thus not appropriate), we used Akaike's
Information Criterium (AIC). Tests for environmental
associations were performed at the amplicon level
(counting the SNPs located within overlapping regions
only once) and including only populations with more than
three gametes sampled (average number of gametes
sampled per population: 7.98 in *P. pinaster* and 9.27 in *P.
halepensis*, see Table S1).

Results

Gene diversity and divergence

Overall primer transfer rates across the two pine species
were high (91.6% for *P. pinaster* and 83.3% for *P. halepensis*)
and we obtained a total of 533 kbp and 380
kbp of sequence in *P. pinaster* and *P. halepensis*,
respectively. All sequences represented putative orthologs
(determined on the basis of similarities with available
sequences in GenBank and construction of phylogenies
including multiple gene members of the target families in
different conifers) except probably for *dhn2-Ps_b* in *P. halepensis*
(see Methods). Maritime pine showed a higher
level of gene diversity than Aleppo pine (Table 2), and this
result held when considering only the nine orthologous
amplicons in the two species (data not shown). In total, 131
SNPs (24 non-synonymous) were detected in *P. pinaster*
(4 to 27 per amplicon) and 65 SNPs (12 non-synonymous)
in *P. halepensis* (1 to 21 per amplicon) (Table 2). There were
only four shared polymorphisms (i.e. 2 % of polymorphic
sites were shared) across species, with two of them (those
within *4cl-Pt*) appearing only in one *P. pinaster* individual.
Number of insertions/deletions was similar between the
two species for the same locus, but these indels were not
shared across species. None of them caused frame shifts.
At the haplotype level, not only the number of haplotypes
substantially differed between the two species (86 for *P. pinaster*
vs. only 59 for *P. halepensis*) but also their
frequency spectrum: *P. pinaster* had multiple common
haplotypes along with some less frequent ones, while *P. halepensis*
generally displayed one major haplotype
together with some low-frequency ones (Figure S2).
Finally, there was also a contrasted pattern of gene
diversity between the regional groups defined by
chloroplast microsatellites in *P. pinaster*: the Western
Mediterranean group displayed higher nucleotide diversity
than the North African group (average per population of
0.0039 vs. 0.0022). The SNP and haplotypes detected
were used to conduct neutrality tests, as well as to detect
associations with environmental variables.

Divergence at the DNA sequence level was examined
using two outgroups: *P. taeda*, a New World pine, and each
of the Mediterranean pines as reference for the other (i.e. *P. halepensis* was used as outgroup for *P. pinaster* and *vice versa*). For both species, nucleotide divergence per site with *P. taeda* (*K*<sub>all</sub> = 0.0397 for *P. pinaster*; *K*<sub>all</sub> = 0.0498 for *P. halepensis*) was higher than with the other Mediterranean species (*K*<sub>all</sub> = 0.0235 for *P. pinaster* vs. *P. halepensis*) except for *dhn2-Pp_b* in *P. pinaster* that displayed a slightly higher divergence between the Mediterranean pines than between each of them with *P. taeda*.

## Neutrality tests

Neutrality tests rejected the null neutral model for distinct genes depending on the species and the geographic group (Table 3). The *D* and *H*<sub>norm</sub> tests identified some genes departing significantly from expectation under the neutral model that are not detected by the more robust *DHEW* test: *lp31-Pt_a* and *lp31-Pt_b* in the North African group for *P. pinaster*, *lp33-Pp* in the Western Mediterranean group for *P. halepensis*. As shown below, these two tests are probably reflecting demographical events and not selective processes. More relevantly, two fragments of the *4cl-Pt* gene were identified as potential targets of selection by *DHEW* tests: *4cl-Pt_b* in the Western Mediterranean group for *P. pinaster* and *4cl-Pt_d* in the North African group for *P. halepensis*. Considering all polymorphic sites or only silent sites did not change the qualitative results of the test (Table 3 and Table S5).

None of the *MK* tests showed any departure from the null hypothesis of an equal ratio of nonsynonymous to synonymous variation within and between species (data not shown). The *ML-HKA* tests based on silent sites (with *P. taeda* as outgroup) revealed uncoupled patterns of polymorphism and divergence in only one amplicon in *P. pinaster* in the Western Mediterranean group (*dhn2-Pp_b*), while in *P. halepensis* one amplicon in this same group (*lp33-Pp*) and one in the North African group (*4cl-Pt_b*) were significant (Table 3). These results still stood when analyses were restricted to only one amplicon per gene (in order to avoid potential linkage between fragments from the same gene) and when fragments potentially under selection were not used as control loci (data not shown). When all substitution sites were included in the analyses, three more genes were significant for *P. halepensis* in the Western Mediterranean group (*dhn2-Ps_a*, *4cl-Pt_a*, *4cl-Pt_b*) (Table S5). In contrast, when Mediterranean pines were used as outgroup for each other, patterns of polymorphism and divergence either considering all sites (Table S5) or silent sites (Table 3) seemed to evolve neutrally and a similar trend of the maximum likelihood estimate of the selection parameter (*ML-HKA's K*) was observed in the two series of analyses.

## Environmental associations

Environmental associations were examined at both the SNP and the haplotype level, and, overall, similar associations were detected. Nonetheless, a few

correlations differed, with some found only at the haplotype level, which highlights how these two levels complement each other as haplotypes may reflect interactions among linked mutations. Within the two species, all significant associations detected between SNPs or haplotypes and climatic variables involved temperature indices as environmental variables.

As many as 23 significant associations were initially found in *P. pinaster* (non-corrected model as provided by SAM, see Methods). Only three associations remained after integrating neutral marker PCs as covariates, two of them with spatial variables (one with altitude, one with latitude) and one (*4cl-Pt_c_229*) with the minimum temperature of the coldest month (MTCM) (Figure 2; Table S6). At the haplotype level, initially 32 associations were found in maritime pine but only six remained in the corrected model: two (*4cl-Pt_c_A* and *4cl-Pt_c_B*) with altitude and one (*dhn5-Ps_C*) with annual mean temperature (AMT), one (*dhn5-Ps_E*) with maximum temperature of the warmest month (MTWM), and two (*4cl-Pt_c_A* and *4cl-Pt_c_B*) with minimum temperature of the coldest month (MTCM) (Figure 2; Table S7). In the case of *P. halepensis*, only one association was found by the corrected models, both at the SNP (between *4cl-Pt_c_131* and temperature seasonality, TS) and at the haplotype (*4cl-Pt_c_B* and TS) level (Figure 2; Tables S6 and S7).

In summary, association analysis identified two loci exhibiting significant correlations with temperature indices: *4cl-Pt_c* that is associated with MTCM in *P. pinaster* and with TS in *P. halepensis* at both the SNP and haplotype levels, and *dhn5-Ps* that is associated with AMT and MTWM at the haplotype level in *P. pinaster* only.

## Discussion

In this study, we assessed the impact of natural selection on the same set of candidate genes related to drought tolerance in two widespread Mediterranean pine species. Our results revealed distinct selection patterns according to species, geographic regions and loci. Below, we discuss these findings in the light of the history of each species and the specificities of each of the methods used to reveal footprints of selection.

### Neutrality tests

Neutrality tests examining selection events based on site- and haplotype-frequency spectrum identified distinct genes potentially targeted by selection. The *D* and/or *H*<sub>norm</sub> tests detected three loci departing from neutrality that were not detected by the *DHEW* compound tests. These results may be explained by the sensitivity of the *D* and *H*<sub>norm</sub> tests to demographic factors and different degrees of background selection (Zeng *et al.* 2006). Especially, Tajima's *D* is sensitive to background selection and population growth while the *H*<sub>norm</sub> is more sensitive to population shrinkage and subdivision. In a previous
study in *P. halepensis* based on ten candidate genes, we have shown that this species has undergone historical bottlenecks and that the western populations of the species harbored some signatures of this demographic event (Grivet *et al.* 2009). In addition, based on the present data set and extensive coalescent simulation, we have found that observed values of *D* and *H*<sub>norm</sub> for *P. pinaster* both in the Western Mediterranean and North African groups reject the standard neutral model and suggest past bottlenecks in this species too (see Figure S3). To test the robustness of the *DHEW* compound test, neutral coalescent simulations under realistic bottlenecks scenarios were further simulated in *P. pinaster* (1,000 coalescent simulations per scenario) and the significant values for this test were recorded. In all cases, the number of significant tests obtained was equal or lower than the number expected by chance (see Figure S3), highlighting the robustness of the *DHEW* test in the presence of bottleneck events and geographic population structure of magnitudes similar to the ones found in maritime and Aleppo pines. In addition to our simulations, the *DHEW* test has been shown to be robust to recombination and to have high sensitivity to detect positive selection (Zeng *et al.* 2007), as it combines powerful (and insensitive to recombination) haplotype-frequency spectrum tests (i.e the *EW* test) with SFS statistics, such as the *D* test that maintains power to detect positive selection across a wide period of time (Zeng *et al.* 2006; Zhai, Nielsen and Slatkin 2009). Only two on-going events of selection were detected by the *DHEW* test, one in *P. pinaster* and one in *P. halepensis*. Both cases involved the *4cl-Pt* gene (albeit in different geographical regions and distinct amplicons), suggesting a potential role of *4cl* in local adaption in both pine species (see below).

Among the tests examining polymorphism within species and divergence among species, the *MK* test did not detect any loci under selection. While this test is somewhat robust to demography, it seems less powerful in detecting positive selection than *HKA* type tests (Zhai, Nielsen, and Slatkin 2009). The *ML-HKA* tests for both silent and all sites, using *P. taeda* as outgroup, identified three amplicons with uncoupled levels of polymorphism and divergence in the two pines. Two genes (*dhn2-Pp_b* in *P. pinaster* in the Western Mediterranean group and *4cl-Pt_b* in *P. halepensis* in the North African group) showed high levels of diversity compared to divergence (selection parameter higher than one), a pattern compatible with balancing selection. One other gene (*lp33-Pp* in the Western Mediterranean group in *P. halepensis*) showed low diversity compared to divergence (selection parameter lower than one), which could reflect the transient reduction in variability occurring during a selective sweep. Three other genes showed uncoupled levels of polymorphism and divergence in *P. halepensis* when examining all sites, but not when examining silent sites alone. Two of these genes (*4cl-Pt_a* and *4cl-Pt_b* in

the Western Mediterranean group) had an excess of nucleotide diversity, which could reflect recent balancing selection acting on non-synonymous sites, while for one gene (*dhn2-Pp_a* in the Western Mediterranean group) purifying selection reducing diversity in non-synonymous sites could explain the observed pattern. Evidence for extensive purifying selection in conifers comes from an approximately four-fold nucleotide diversity at synonymous compared to non-synonymous sites in most species (see Table 1 in González-Martínez *et al.* 2010), including maritime (dN/dS = 0.169) and Aleppo (dN/dS = 0.344) pines (this study). Although demographical events can result in a significant *ML-HKA* test, its multilocus nature combined with the *HKA* framework should produce a more robust test than those based on comparing different aspects of polymorphism at a single locus (Wright and Charlesworth 2004).

In order to gain power by increasing sample size while taking the specific history of each species into account, we performed the neutrality tests in groups of populations that present similar evolutionary history (i.e. North African and Western Mediterranean groups). However, grouping populations may bias the output of the SFS-based neutrality tests as it can increase the proportion of singletons, leading to an excess of low frequency variants and thus to more negative *D* (Städler *et al.* 2009). In our study, populations were grouped according to homogeneous gene pools (see Bucci *et al.* 2007 for *P. pinaster* and Grivet *et al.* 2009 for *P. halepensis*), the consequence of which should be that no substantial differences in the number of rare variants at the population and group levels are found. The extensive gene flow normally found in conifers with large and continuous distribution would also support this approach. Nevertheless, to check that our results were not biased because of population grouping, we performed neutrality tests at the population level too and found that values of *D* test at the group level were not more negative than at the single population level (data not shown). Moreover, results of the *H*<sub>norm</sub> and *DHEW* tests should not be affected by the accumulation of singletons as the first is based on intermediary and high frequency variants while the second integrates all variants of the site frequency spectrum. Thus, our analyses appear reasonably robust to population grouping.

Selective events appear to have affected distinct genes in *P. pinaster* and in *P. halepensis* despite their overlapping environment and close phylogenetic relationship, a result that may be connected to the different histories of the two pines that have likely resulted in different selective pressures. In *P. pinaster*, selective events were detected only within the Western Mediterranean group, while in *P. halepensis* footprints of selection were identified both in the Western Mediterranean and North African groups. These results point to different geographical selection pressures that may have led to the process of regional
adaptation of the pines (see Barbéro *et al.* 1998 for a description of the different Mediterranean environmental zones; and Gómez and Zamora 2000; Nakazato, Bogonovich, and Moyle 2008; Montesinos *et al.* 2009 for some examples on adaptive variation across heterogeneous environments). It is also noticeable that *P. pinaster* had lower levels of nucleotide variation within the North African range. This fact can be attributed to population history and/or interpreted in the light of more extreme environmental conditions constraining population sizes and/or resulting in stronger selection in this range.

## Environmental associations

Environmental associations identified two loci that were correlated with temperature, suggesting the importance of this climatic variable as selective agent (Saxe *et al.* 2001; Jump *et al.* 2006). One of these genes was common for both pines (4c1-Pt correlated with MTCM in *P. pinaster* and with TS in *P. halepensis*), while the other gene association was only significant in *P. pinaster* (*dhn5-Ps* correlated with AMT and MTWM). Some of the correlations between genetic and climatic data in *P. pinaster* were due basically to the extreme values of a few North African populations (Sidi Meskour and Tamrabta for 4c1-Pt_c; Tabarka for *dhn5-Ps_C*; see Figures 2 and 3) and may not represent true adaptive responses to environmental gradients but local adaptation to particular environments or genetic drift due to population isolation (see references in Alleaume-Benharira, Pen, and Ronce 2006; Rosenblum, Hickerson, and Moritz 2007). In contrast, the other associations (4c1-Pt_c in *P. halepensis* and *dhn5-Ps_E* in *P. pinaster*) showed more robust patterns. For instance, *dhn5-Ps_E* tended to be absent from populations characterized by the highest MTWM (Arenas de San Pedro in central Spain, Tamrabta in Moroccan Middle Atlas and Tabarka in coastal Tunisia) regardless of their geographical location (Table 1 and Figure 3).

Environmental correlations with allelic variation can be spuriously inflated by neutral processes that may also generate genetic clines, such as population history or population genetic structure. This is particularly true for forest trees that may have followed a postglacial colonization pathway overlapping temperature and rainfall clines. Here, to control for confounding associations between candidate genes and environmental data, we included variation of chloroplast microsatellites as a covariate in the multivariate logistic regressions. In *P. pinaster*, the use of cpSSRs has allowed a more accurate description of genetic structure (Bucci *et al.* 2007) than previous studies based on biochemical markers, such as terpenes or allozymes (e.g. Baradat and Marpeau-Bezard 1988; Salvador *et al.* 2000). Recently, nuclear molecular data for a geographically limited set of populations (nuSSRs and SNPs) (Eveno *et al.* 2008; our unpublished results) have identified similar gene pools as Bucci *et al.* (2007). More extensive sampling with highly-polymorphic

biparentally-inherited nuclear markers could, in principle, reveal a more complex spatial genetic pattern affecting the outcome of some of the multivariate logistic regressions reported in this study. However, because the available chloroplast data set we used strongly reflects the species' history, they are expected to 'correct' for the presence of overall neutral genetic gradients. Accordingly, including these markers as covariates identified a substantial fraction of the correlations initially retained as significant as false positives (83.6% in *P. pinaster* and 85.7% in *P. halepensis*) (Tables S6 and S7).

## The dehydrin gene family

Neutrality tests and environmental associations both point to dehydrins as potential targets of natural selection in *P. pinaster*. Some dehydrins were also suggested to be under selection in a study based on detection of outlier loci in this species (Eveno *et al.* 2008). Dehydrins displayed also non-neutral patterns of nucleotide diversity in *Pinus sylvestris* populations showing divergence for cold tolerance in Europe (Wachowiak *et al.* 2009) and were associated with carbon isotope discrimination (and, thus, potentially with drought tolerance) in *P. taeda* (González-Martínez *et al.* 2008). Altogether these studies suggest the involvement of dehydrins in the adaptive response of pines to abiotic stress.

Dehydrins are part of a relatively small multigene family of intracellular stabilizers that plays a major role in cell protection against desiccation. These proteins are produced in response to any type of stress that causes dehydration at the cellular level, such as cold, drought or salinity (Close 1997). Changes in dehydrin gene expression have been reported in response to drought and/or cold stress in many plants such as cowpea (Ismail, Hall, and Close 1999), barley (Suprunova *et al.* 2004), wheat (Lopez *et al.* 2001), apple (Wisniewski *et al.* 2008), tomato (Weiss and Egea-Cortines 2009) and blueberry (Panta, Rieger, and Rowland 2001); while transgenic experiments confirmed the role of dehydrin genes in enhancing tolerance to drought or freezing stress in plants (e.g. in Arabidopsis, Puhakainen *et al.* 2004). In conifers, at least eight dehydrin genes have been identified in *Pinus sylvestris* (Joosen *et al.* 2006) and *Picea abies* (Yakovlev *et al.* 2008), and dehydrin expression has been shown to increase under wounding, cold and drought stress (Richard *et al.* 2000; Watkinson *et al.* 2003). Thus, it is not surprising to find in our study significant associations between dehydrins and temperatures, as critical low temperatures can cause tissue injury while high temperatures accompany dehydration, both stimulating the accumulation of dehydrins (Lewitt 1980; Ingram and Bartels 1996; Rizhsky, Liang, and Mittler 2002). Although our results, and others in *P. pinaster*, suggest that selection may have acted on some of the dehydrins tested, further genetic association and functional studies are necessary to confirm the role of dehydrins in local adaptation of this Mediterranean pine.
## The 4-coumarate: CoA ligase (4cl) gene

There was converging evidence both from neutrality tests (*DHEW* in *P. pinaster* and both *DHEW* and *ML-HKA* in *P. halepensis*) and environmental associations that the 4-coumarate: CoA ligase (4cl) gene may be under selection in the two Mediterranean pines studied. Preliminary results on phenotypical association corrected by neutral gradients in *P. pinaster* also point at 4cl-Pt. Indeed, multiple regressions reveal significant associations between polymorphism in 4cl-Pt_c at the SNP and haplotype levels and total height in populations growing in sites characterized by dry and intermediate humidity (Table S8). The 4cl family has been extensively studied in plants, where it is encoded by four to five genes in the fully sequenced genomes of Arabidopsis, rice and poplar (reviewed in Souza *et al.* 2008). The 4cl gene is involved, among other processes, in the production of basic enzymes of the phenylpropanoid metabolism that are important metabolites acting as protectants against biotic and abiotic stresses (Rani *et al.* 2009). The 4cl gene also encodes key enzymes in the biosynthesis of lignin and several studies have demonstrated its involvement in plant growth (e.g. Yun *et al.* 2005; Wagner *et al.* 2009; Yun *et al.* 2009). Thus, we expect that changes in 4cl function would have significant repercussions on tree physiology and morphology. Implication of 4cl in pine morphology and physiology has been shown in gene association studies (for *P. taeda* see González-Martínez *et al.* 2007b) as well as in gene suppression studies (for *P. radiata* see Wagner *et al.* 2009).

## Neutrality tests versus Environmental correlations

The different approaches used in this study suggested different loci under selection within each of the two Mediterranean pine species. This fact has to be connected to the specificities of the statistical methods used, whose performance has been recently studied under various demographic and recombination scenarios (Zeng *et al.* 2007; Zeng, Shi, and Wu 2007; Ramírez-Soriano *et al.* 2008; Zhai, Nielsen, and Slatkin 2009). Within-species site and haplotype frequency spectrum-based methods are suitable to detect ongoing or recently fixed selective sweeps - *EW* tending to be more powerful around the time when a selected mutation reaches fixation (Zeng, Shi, and Wu 2007). In contrast, tests comparing levels of (within-species) polymorphism and (among-species) divergence are able to detect the cumulative effects of positive selection events over a wider evolutionary scale. Within this category, the *ML-HKA* test is suitable to quantify the amount of selection in the genome and, due to its multilocus nature, it is expected to be more robust to changes in population size such as population bottlenecks and expansion than traditional tests (Wright and Charlesworth 2004). Simulations have shown that the *HKA* test is relatively insensitive to change of

divergence time as most of its power comes from the transient reduction in variability occurring during selective sweeps (Zhai, Nielsen, and Slatkin 2009). However, we only obtained significant tests when using the less closely-related outgroup (i.e. a New World pine, *P. taeda*). This points to either lower power when divergence time from the outgroup is low or, alternatively, to selection events that took place before the split of the two Mediterranean pines considered. In our study, none of the *EW* neutrality tests were significant (data not shown), while both the *DHEW* and the *ML-HKA* tests revealed distinct genes potentially under selection. As a consequence, potential selective events detected by the neutrality tests assessed in this study would correspond to either on-going or relatively old selective events (i.e. before the diversification of the Mediterranean pines in the Miocene, ~10 Ma; Gernandt *et al.* 2008).

Last generation tests, such as *DHEW* or *ML-HKA*, are powerful for detecting positive selection and are relatively insensitive to other evolutionary forces, but still do not integrate recombination rate, a factor that can produce substantial biases (Nielsen *et al.* 2007; Ramírez-Soriano *et al.* 2008), nor provide insights on selection drivers. Methods based on correlation between genetic and environmental data are appealing because they aim at understanding species-specific adaptations and processes that connect an organism to its environment by looking at allelic frequencies or the genetic structure of populations (Foll and Gaggiotti 2006; Joost, Kalbermatten, and Bonin 2008). These correlation methods are also more appealing than methods based on detection of outliers (e.g. Eveno *et al.* 2008) as they target particular selection drivers and provide a hypothesis testing framework. In addition, they are well suited for high-throughput genotyping even for non-model species (Namroud *et al.* 2008; Eckert *et al.* forthcoming). However, correlation methods are not without limitations, as it may be challenging to find the environmental factors that are relevant for each species adaptation and associations have to be controlled for historical and demographical processes, in particular for allelic clines produced by postglacial migrations in Europe and the Americas.

Together, neutrality tests and environmental association approaches complement each other by looking at different evolutionary scales and types of selection. In our study, they detected a relatively-high number of genes showing non-neutral patterns of evolution, a result that can be attributed to a selection of candidate genes based on earlier studies (González-Martínez *et al.* 2006; Eveno *et al.* 2008; Grivet *et al.* 2009). This fact supports a candidate gene approach with targeted genes, at least for organisms that have large genomes (e.g. conifers), which, so far, preclude dense genome-wide sampling.
## Conclusion

Our pluralistic approach revealed a dynamic action of natural selection in space and time within maritime and Aleppo pines (Felsenstein 1976; see Vasemagi and Primmer 2005 for other examples). Selection events along with environmental associations were detected; some of these events differ between the two species reflecting individual histories (recolonization, demography, adaptation) while others are shared, which translates partly as a common history of these closely-related and partially sympatric Mediterranean pines. The simultaneous search for patterns of selection in two closely-related species can help understanding the evolutionary forces responsible for adaptive responses, and thus provides an effective way of assessing the degree of local adaptation, a key factor to integrate in future management and conservation strategies (Wright and Gaut 2005).

## Acknowledgements

We are grateful to C. García-Barriga for her assistance in producing sequences of *P. pinaster* candidate genes, to C. Collada, M.T. Cervera, M.A. Guevara, G. Gill, G.R. Brown and D.B. Neale for providing some of the primers, as well as to A.I. de-Lucas and C. Ordoñez for supplying some of the *P. pinaster* seeds. We thank J. Gonzalo for kindly providing the climatic data from his functional phytoclimatic model, as well as S.E. Ramos-Onsins for his scripts to compute the diversity parameters and K. Zeng for his scripts to compute the *DHEW* neutrality test. Thanks are extended to M. Robson, M.A. Zavala and J. Climent for helpful discussion on the study. We are also thankful to A.J. Eckert for thoughtful comments on the manuscript, N. Takebayashi and two anonymous reviewers for constructive comments and suggestions, and P.C. Grant for revising the English. The EUFORGEN program (Bioversity International) provided the distribution maps used in Fig. 1.

## Funding

This work was supported by the EU EVOLTREE Network of Excellence, the Collaborative Project on 'Conservation of Forest Genetic Resources' between the Spanish Ministry of Environment and INIA (AEG06-054) and Projects CGL2008-05289-C02-02 /BOS (VaMPiro) and EUI2008-03713 (BiodivERsA-LinkTree) from the Spanish Ministry of Science and Innovation, the later through the FP6 BiodivERsA Eranet. TB acknowledges support from the Danish Research Council (FNU) through a Steno fellowship.

## LITERATURE CITED

Alleaume-Benharira M, Pen IR, Ronce O. 2006. Geographical patterns of adaptation within a species' range: interactions between drift and gene flow. J Evol Biol. 19:203-215.

Aranda I, Alía R, Ortega U, Dantas AK, Majada J. 2009. Intra-specific variability in biomass partitioning and carbon isotopic discrimination under moderate drought stress in seedlings from *Pinus pinaster* populations. Tree Genet Genomes 6:169-178.

Atzmon N, Moshe Y, Schiller G. 2004. Ecophysiological response to severe drought in *Pinus halepensis* Mill. trees of two provenances. Plant Ecol. 171:15-22.

Baradat PH, Marpeau-Bezard A. 1988. Le pin maritime, *Pinus pinaster* Ait. Biologie et génétique des terpènes pour la connaissance et l'amélioration de l'espèce. PhD Thesis, University of Bordeaux I, France.

Barbéro M, Loisel R, Quézel P, Richardson DM, Romane F. 1998. Pines of the Mediterranean Basin. In Richardson DM, editor. Ecology and Biogeography of *Pinus*. Cambridge University Press, Cambridge. Pp. 153-170.

Biswas S, Akey JM. 2006. Genomic insights into positive selection. Trends Genet. 22:437-446.

Bucci G, Anzidei M, Madaghiele A, Vendramin GG. 1998. Detection of haplotypic variation and natural hybridization in halepensis-complex pine species using chloroplast simple sequence repeat (SSR) markers. Mol Ecol. 7:1633-1643.

Bucci G, González-Martínez SC, Le Provost G, Plomion C, Ribeiro MM, Sebastiani F, Alía R, Vendramin GG. 2007. Range-wide phylogeography and gene zones in *Pinus pinaster* Ait. revealed by chloroplast microsatellite markers. Mol Ecol. 16:2137-2153.

Burban C, Petit RJ. 2003. Phylogeography of maritime pine inferred with organelle markers having contrasted inheritance. Mol Ecol. 12:1487-1495.

Chambel MR, Climent J, Alía R. 2007. Divergence among species and populations of Mediterranean pines in biomass allocation of seedlings grown under two watering regimes. Ann For Sci. 64:87-97.

Climent J, Prada MA, Calama R, Chambel MR, De Ron DS, Alía R. 2008. To grow or to seed: Ecotypic variation in reproductive allocation and cone production by young female Aleppo pine (*Pinus halepensis*, Pinaceae). Am J Bot. 95:833-842.

Close T. 1997. Dehydrins: a commonality in the response of plants to dehydration and low temperature. Physiol Plant. 100:291-296.

Coop G, Pickrell JK, Novembre J, Kudaravalli S, Li J, Absher D, Myers RM, Cavalli-Sforza LL, Feldman MW, Pritchard JK. 2009. The role of geography in human adaptation. PLoS Genetics 5(6): e1000500.

Costa P, Bahrman N, Frigerio JM, Kremer A, Plomion C. 1998. Water-deficit-responsive proteins in maritime pine. Plant Mol Biol. 38:587-596.

Dellaporta SL, Wood J, Hicks JB. 1983. A plant DNA minipreparation: version II. Plant Mol Biol Rep. 1:19-21.

Eckert AJ, van Heerwaarden J, Wegrzyn JL, Nelson CD, Ross-Ibarra J, González-Martínez SC, Neale DB. Forthcoming. Patterns of population structure and environmental associations to aridity across the range of loblolly pine (*Pinus taeda* L., Pinaceae). Genetics.

Eveno E, Collada C, Guevara MA, et al. (11 co-authors). 2008. Contrasting patterns of selection at *Pinus pinaster* Ait. drought stress candidate genes as revealed by genetic differentiation analyses. Mol Biol Evol. 25:417-437.
Felsenstein J. 1976. Theoretical population-genetics of variable selection and migration. Annu Rev Genet.10:253-280.

Foll M, Gaggiotti O. 2006. Identifying the environmental factors that determine the genetic structure of populations. Genetics 174:875-891.

Frankel N, Carrari F, Hasson E, Iusem ND. 2006. Evolutionary history of the Asr gene family. Gene 378:74-83.

Gernandt DS, Lopez GG, Garcia SO, Liston A. 2005. Phylogeny and classification of Pinus. Taxon 54:29-42.

Gernandt DS, Magallon S, Lopez GG, Flores OZ, Willyard A, Liston A. 2008. Use of simultaneous analyses to guide fossil-based calibrations of *Pinaceae* phylogeny. Int J Plant Sci.169:1086-1099.

Giorgi F, Lionello P. 2008. Climate change projections for the Mediterranean region. Glob Planet Change 63:90-104.

Gómez JM, Zamora R. 2000. Spatial variation in the selective scenarios of *Hormathophylla spinosa* (Cruciferae). Am Nat. 155:657-668.

González-Martínez SC, Dillon S, Garnier-Géré P, et al. (16 co-authors). Forthcoming. Patterns of nucleotide diversity and association mapping In Kole C, editor. Genomics of Conifers. Science Publishers, Inc., New Hampshire.

González-Martínez SC, Ersoz E, Brown GR, Wheeler NC, Neale DB. 2006. DNA sequence variation and selection of tag single-nucleotide polymorphisms at candidate genes for drought-stress response in *Pinus taeda* L. Genetics 172:1915-1926.

González-Martínez SC, Gómez A, Carrión JS, Agúndez D, Alía R, Gil L. 2007a. Spatial genetic structure of an explicit glacial refugium of maritime pine (*Pinus pinaster* Aiton) in southeastern Spain. In Weiss S, Ferrand N, editors. Phylogeography in southern European refugia: evolutionary perspectives on the origins of European biodiversity, London. Pp. 257-269.

González-Martínez SC, Huber D, Ersoz E, Davis JM, Neale DB. 2008. Association genetics in *Pinus taeda* L. II. Carbon isotope discrimination. Heredity 101:19-26.

González-Martínez SC, Wheeler NC, Ersoz E, Nelson CD, Neale DB. 2007b. Association genetics in *Pinus taeda* L. I. Wood property traits. Genetics 175:399-409.

Gonzalo J. 2007. Phytoclimatic analysis of the Spanish Peninsula. Update and geostatistical analysis. PhD Thesis, University of Valladolid, Spain.

Gram WK, Sork VL. 2001. Association between environmental and genetic heterogeneity in forest tree populations. Ecology 82:2012-2021.

Grivet D, Sebastiani F, González-Martínez SC, Vendramin GG. 2009. Patterns of polymorphism resulting from long-range colonization in the Mediterranean conifer Aleppo pine. New Phytol. 184:1016-1028.

Hancock AM, Witonsky DB, Gordon AS, Eshel G, Pritchard JK, Coop G, Di Rienzo A. 2008. Adaptations to climate in candidate genes for common metabolic disorders. PLoS Genetics 4(2): e32.

Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 25:1965-1978.

Hudson RR, Kreitman M, Aguade M. 1987. A test of neutral molecular evolution based on nucleotide data. Genetics 116:153-159.

Ingram J, Bartels D. 1996. The molecular basis of dehydration tolerance in plants. Annu Rev Plant Physiol Plant Mol Biol. 47:377-403.

Ismail AM, Hall AE, Close TJ. 1999. Allelic variation of a dehydrin gene cosegregates with chilling tolerance during seedling emergence. Proc Natl Acad Sci USA. 96:13566-13570.

Joosen RVL, Lammers M, Balk PA, Bronnum P, Konings MCJM, Perks M, Stattin E, Van Wordragen MF, van der Geest AHM. 2006. Correlating gene expression to physiological parameters and environmental conditions during cold acclimation of *Pinus sylvestris*, identification of molecular markers using cDNA microarrays. Tree Physiol. 26:1297-1313.

Joost S, Kalbermatten M, Bonin A. 2008. Spatial analysis method (SAM): a software tool combining molecular and environmental data to identify candidate loci for selection. Mol Ecol Resources 8:957-960.

Jump AS, Hunt JM, Martínez-Izquierdo JA, Peñuelas J. 2006. Natural selection and climate change: temperature-linked spatial and temporal trends in gene frequency in *Fagus sylvatica*. Mol Ecol. 15:3469-3480.

Keller SR, Sowell DR, Neiman M, Wolfe LM, Taylor DR. 2009. Adaptation and colonization history affect the evolution of clines in two introduced species. New Phytol. 183:678-690.

Kimura M, Ohta T. 1974. Some principles governing molecular evolution. Proc Natl Acad Sci USA. 71:2848-2852.

Lewitt J. 1980. Responses of plants to environmental stresses, Vol I. Chilling, freezing and high temperature stresses. Academic Press, New York, USA.

Librado P, Rozas J. 2009. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25:1451-1452.

Lopez CG, Banowetz G, Peterson CJ, Kronstad WE. 2001. Differential accumulation of a 24-kd dehydrin protein in wheat seedlings correlates with drought stress tolerance at grain filling. Hereditas 135:175-181.

Manel S, Schwartz MK, Luikart G, Taberlet P. 2003. Landscape genetics: combining landscape ecology and population genetics. Trends Ecol Evol. 18:189-197.

McDonald JH, Kreitman M. 1991. Adaptive protein evolution at the ADH locus in *Drosophila* Nature 351:652-654.

Montesinos A, Tonsor SJ, Alonso-Blanco C, Picó FX. 2009. Demographic and genetic patterns of variation among populations of *Arabidopsis thaliana* from contrasting native environments. PLoS One 4:e7213.

Morgante M, Felice N, Vendramin GG. 1998. Analysis of hypervariable chloroplast microsatellites in *Pinus halepensis* reveals a dramatic genetic bottleneck. In Karp A, Isaac PG, Ingram DS, editors. Molecular tools for screening biodiversity. Chapman & Hall, London, UK. Pp. 407-412.

Nakazato T, Bogonovich M, Moyle LC. 2008. Environmental factors predict adaptive phenotypic differentiation within and between two wild Andean tomatoes. Evolution 62:774-792.

Namroud MC, Beaulieu J, Juge N, Laroche J, Bousquet J. 2008. Scanning the genome for gene single nucleotide polymorphisms involved in adaptive population differentiation in white spruce. Mol Ecol. 17:3599-3613.
Nei M. 1987. Molecular evolutionary genetics. Columbia Univ. Press, New York.

Nguyen A, Lamant A. 1989. Variation in growth and osmotic regulation of roots of water-stressed maritime pine (Pinus pinaster Ait) provenances. Tree Physiol. 5:123-133.

Nielsen R. 2005. Molecular signatures of natural selection. Annu Rev Genet. 39:197-218.

Nielsen R, Hellmann I, Hubisz M, Bustamante C, Clark AG. 2007. Recent and ongoing selection in the human genome. Nat Rev Genet. 8:857-868.

Panta GR, Rieger MW, Rowland LJ. 2001. Effect of cold and drought stress on blueberry dehydrin accumulation. J Horticult Sci Biotechnol. 76:549-556.

Parisod C, Christin PA. 2008. Genome-wide association to fine-scale ecological heterogeneity within a continuous population of *Biscutella laevigata* (Brassicaceae). New Phytol. 178:436-447.

Petit RJ, Hampe A, Cheddadi R. 2005. Climate changes and tree phylogeography in the Mediterranean. Taxon 54:877-885.

Puhakainen T, Hess MW, Mäkelä P, Svensson J, Heino P, Palva ET. 2004. Overexpression of multiple dehydrin genes enhances tolerance to freezing stress in *Arabidopsis*. Plant Mol Biol. 54:743-753.

Quézel P. 2000. Taxonomy and biogeography of Mediterranean pines. In Ne'eman G, Trabaud L, editors. Ecology, biogeography and management of *Pinus halepensis* and *P. brutia* forest ecosystems in the Mediterranean Bassin. Backhuys Publishers, Leiden, The Netherlands. Pp. 1-12.

Ramírez-Soriano A, Ramos-Onsins SE, Rozas J, Calafell F, Navarro A. 2008. Statistical power analysis of neutrality tests under demographic expansions, contractions and bottlenecks with recombination. Genetics 179:555-567.

Rani A, Singh K, Sood P, Kumar S, Ahuja PS. 2009. p-Coumarate:CoA ligase as a key gene in the yield of catechins in tea [*Camellia sinensis* (L.) O. Kuntze]. Funct Integr Genomics 9: 271-275.

Resch T. 1974. Essai de distinction des races majeurs de *Pinus pinaster*. Annales de la Recherche Forestière au Maroc 14:91-102.

Richard S, Morency MJ, Drevet C, Jouanin L, Séguin A. 2000. Isolation and characterization of a dehydrin gene from white spruce induced upon wounding, drought and cold stresses. Plant Mol Biol. 43:1-10.

Richardson DM, Bond WJ. 1991. Determinants of plant-distribution - Evidence from pine invasions. Am Nat. 137:639-668.

Rizhsky L, Liang HJ, Mittler R. 2002. The combined effect of drought stress and heat shock on gene expression in tobacco. Plant Physiol. 130:1143-1151.

Rosenblum EB, Hickerson MJ, Moritz C. 2007. A multilocus perspective on colonization accompanied by selection and gene flow. Evolution 61:2971-2985.

Salvador L, Alía R, Agúndez D, Gil L. 2000. Genetic variation and migration pathways of maritime pine (*Pinus pinaster* Ait) in the Iberian peninsula. Theor Appl Genet. 100:89-95.

Sathyan P, Newton RJ, Loopstra CA. 2005. Genes induced by WDS are differentially expressed in two populations of

Aleppo pine (*Pinus halepensis*). Tree Genet Genomes 1:166-173.

Savolainen O, Pyhäjärvi T, Knürr T. 2007. Gene flow and local adaptation in trees. Annu Rev Ecol Evol Syst. 38:595-619.

Saxe H, Cannell MGR, Johnsen B, Ryan MG, Vourlitis G. 2001. Tree and forest functioning in response to global warming. New Phytol. 149:369-399.

Soto A, Robledo-Arnuncio JJ, González-Martínez SC, Smouse PE, Alía R. 2010. Climatic niche and neutral genetic diversity of the six Iberian pine species: a retrospective and prospective view. Mol Ecol. 19:1396-1409.

Souza CD, Barbazuk B, Ralph SG, Bohlmann J, Hamberger B, Douglas CJ. 2008. Genome-wide analysis of a land plant-specific acyl: coenzymeA synthetase (ACS) gene family in *Arabidopsis*, poplar, rice and *Physcomitrella*. New Phytol. 179:987-1003.

Städler T, Haubold B, Merino C, Stephan W, Pfaffelhuber P. 2009. The impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided populations. Genetics 182:205-216.

Storz JF. 2002. Contrasting patterns of divergence in quantitative traits and neutral DNA markers: analysis of clinal variation. Mol Ecol. 11:2537-2551.

Suprunova T, Krugman T, Fahima T, Chen G, Shams I, Korol A, Nevo E. 2004. Differential expression of dehydrin genes in wild barley, *Hordeum spontaneum*, associated with resistance to water deficit. Plant Cell Environ. 27:1297-1308.

Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585-595.

Tapias R, Climent J, Pardos JA, Gil L. 2004. Life histories of Mediterranean pines. Plant Ecol. 171:53-68.

Vasemagi A, Primmer CR. 2005. Challenges for identifying functionally important genetic variation: the promise of combining complementary research strategies. Mol Ecol. 14:3623-3642.

Voltas J, Chambel M, Prada M, Ferrio J. 2008. Climate-related variability in carbon and oxygen stable isotopes among populations of Aleppo pine grown in common-garden tests. Trees-Struct Funct. 22:759-769.

Wachowiak W, Balk PA, Savolainen O. 2009. Search for nucleotide diversity patterns of local adaptation in dehydrins and other cold-related candidate genes in Scots pine (*Pinus sylvestris* L.). Tree Genet Genomes 5:117-132.

Wagner A, Donaldson L, Kim H, Phillips L, Flint H, Steward D, Torr K, Koch G, Schmitt U, Ralph J. 2009. Suppression of 4-Coumarate-CoA ligase in the coniferous gymnosperm *Pinus radiata*. Plant Physiol. 149:370-383.

Watkinson JI, Sioson AA, Vasquez-Robinet C, et al. (14 co-authors). 2003. Photosynthetic acclimation is reflected in specific patterns of gene expression in drought-stressed loblolly pine. Plant Physiol. 133: 1702-1716.

Watterson GA. 1975. On the number of segregating sites in genetical models without recombination. Theor Pop Biol. 7: 256-276.

Watterson GA. 1978. The homozygosity test of neutrality. Genetics 88:405-417.

Weiss J, Egea-Cortines M. 2009. Transcriptomic analysis of cold response in tomato fruits identifies dehydrin as a marker of cold stress. J Appl Genetics 50:311-319.
Wisniewski M, Bassett C, Norelli J, Macarisin D, Artlip T, Gasic K, Korban S. 2008. Expressed sequence tag analysis of the response of apple (Malus x domestica 'Royal Gala') to low temperature and water deficit. Physiol Plant. 133:298-317.

Wright SI, Charlesworth B. 2004. The HKA test revisited: a Maximum-Likelihood-Ratio test of the standard neutral model. Genetics 168:1071-1076.

Wright SI, Gaut BS. 2005. Molecular population genetics and the search for adaptive evolution in plants. Mol Biol Evol. 22:506-519.

Yakovlev IA, Asante DKA, Fossdal CG, Partanen J, Junttila O, Johnsen O. 2008. Dehydrins expression related to timing of bud burst in Norway spruce. Planta 228:459-472.

Yun MS, Chen WJ, Deng F, Yogo Y. 2005. Differential properties of 4-coumarate: CoA ligase related to growth suppression by chalcone in maize and rice. Plant Growth Regul. 46:169-176.

Yun MS, Chen WJ, Deng F, Yogo Y. 2009. Selective growth suppression of five annual plant species by chalcone and naringenin correlates with the total amount of 4-coumarate: coenzyme A ligase. Weed Biol Manag. 9:27-37.

Zeng K, Fu Y-X, Shi S, Wu C-I. 2006. Statistical tests for detecting positive selection by utilizing high-frequency variants. Genetics 174:1431-1439.

Zeng K, Mano SH, Shi SH, Wu CI. 2007. Comparisons of site- and haplotype-frequency methods for detecting positive selection. Mol Biol Evol. 24:1562-1574.

Zeng K, Shi S, Wu CI. 2007. Compound tests for the detection of hitchhiking under positive selection. Mol Biol Evol. 24:1898-1908.

Zhai WW, Nielsen R, Slatkin M. 2009. An investigation of the statistical power of neutrality tests based on comparative and population genetic data. Mol Biol Evol. 26:273-283.
FELIX GUGERLI<fnref n="1" />, ROLAND BRANDL<fnref n="2" />, BASTIEN CASTAGNEYROL<fnref n="3" />, ALAIN FRANC<fnref n="3" />, HERVÉ JACTEL<fnref n="3" />, HANS-PETER KOELEWIJN<fnref n="4" />, FRANCIS MARTIN<fnref n="5" />, MARTINA PETER<fnref n="1" />, KARIN PRITSCH<fnref n="6" />, HILKE SCHRÖDER<fnref n="7" />, MARINUS J.M. SMULDERS<fnref n="8" />, ANTOINE KREMER<fnref n="3" />, BIRGIT ZIEGENHAGEN<fnref n="2" /> and EVOLTREE JERA<fnref n="3" /> CONTRIBUTORS*

**Corresponding Author:**

Felix Gugerli
WSL Swiss Federal Research Institute, 8903 Birmensdorf, Switzerland
E-mail: gugerli@wsl.ch)

Understanding the interactions of co-occurring species within and across trophic levels provides key information needed for understanding the ecological and evolutionary processes that underlie biological diversity. As genetics has only recently been integrated into the study of community-level interactions, the time is right for a critical evaluation of potential new, gene-based approaches to studying communities. Next generation molecular techniques, used in parallel with field-based observations and manipulative experiments across spatio-temporal gradients, are key to expanding our understanding of community-level processes. Here, we introduce a variety of "-omics" tools, with recent studies of plant-insect herbivores and of ectomycorrhizal systems providing detailed examples of how next generation approaches can revolutionize our understanding of interspecific interactions. We suggest ways that novel technologies may convert community genetics from a field that relies on correlative inference to one that reveals causal mechanisms of genetic co-variation and adaptations within communities.

**Key words:** ectomycorrhizal symbiosis, forest ecosystems, gene-to-gene interactions, herbivorous insects, population genomics, quantitative trait analysis.

Published in *Molecular Ecology* (2013) 22: 3198-3207, with doi:10.1111/mec.12300

Community genetics aims to understand how genetic variation within and among populations of host species affects the composition of associated organisms interacting with the host (Agrawal 2003; Whitham et al. 2006; Johnson & Stinchcombe 2007; Rowntree et al. 2011; Wymore et al. 2011). Empirical community genetics has been stimulated by pioneering work on poplars (*Populus* spp.), their genotype-based phenotypic variation, and associated communities (Whitham et al. 2006). However, community genetics has hitherto largely remained phenomenological, and the underlying genetic basis and processes involved in the interactions between host and associated organisms have not been studied in detail yet. Given the rapid development of molecular techniques (Rokas & Abbot 2009), it will soon be feasible to characterize the genomes of numerous members of a community. With whole-genome sequences or other types of -omics data at hand (Nadeau & Jiggins 2010), community genetics will be able to establish a solid genetic framework in which to understand the interplay between ecological and evolutionary processes (Rokas &

Abbot 2009). Here, we sketch possible avenues along which research in community genetics may proceed, focussing in particular on how -omics may improve our understanding of the role of gene variants in species interactions. First, we argue for exploring spatio-temporal variation to investigate the fundamental ecological and evolutionary aspects of community genetics. Second, we describe how genomic, transcriptomic, proteomic, and metabolomic research can improve understanding of the interactions between trees as focal species and ectomycorrhizal fungi or herbivorous insects, the key players in forest ecosystems.

## Community genetics in a spatio-temporal perspective

Let us consider populations of a focal species that start to diverge genetically. Genetic drift and/or selection may induce shifts in allele frequencies, leading to changes in the phenotypic traits mediating interactions with associated species that use the focal species as a host.

<footnote n="1">¹ WSL Swiss Federal Research Institute, 8903 Birmensdorf, Switzerland</footnote>

<footnote n="2">² Fachbereich Biologie, Philipps-Universität Marburg, 35032 Marburg, Germany</footnote>

<footnote n="3">³ UMR1202 Biodiversity, Genes & Communities, INRA Pierroton, 33612 Cestas Cedex, France</footnote>

<footnote n="4">⁴ ALTERRA Centre for Ecosystem Studies, Wageningen UR, 6700 AA Wageningen, The Netherlands</footnote>

<footnote n="5">⁵ UMR "Interactions Arbres/Micro-Organismes", INRA Nancy, 54280 Champenoux, France</footnote>

<footnote n="6">⁶ Institute of Soil Ecology, Helmholtz Zentrum München, 85764 Neuherberg, Germany</footnote>

<footnote n="7">⁷ Institute for Forest Genetics, Johann Heinrich von Thuenen-Institute, 22927 Grosshansdorf, Germany</footnote>

<footnote n="8">⁸ Plant Research International, Wageningen UR, 6700 AA Wageningen, The Netherlands</footnote>

<footnote>* The EC-supported Network of Excellence Evoltree (http://www.evoltree.eu) formed a group of scientists involved in and actively contributing to Jointly Executed Research Activities on community genetics in forest ecosystems. In addition to the main authors, the group includes S. Augustin, M. Brändle, C. Burban, J. Burczyk, S. Cavers, I. Chybicki, C. Conord, E. Cremer, J. DeWoody, K. Donges, B. Fady, L. Karlinski, C. Kerdelhué, B. Kieliszewska-Rokicka, G. Kost, M. Kulczyk-Skrzeszewska, F. Lakatos, F. Lefèvre, S. Liepelt, S. Oddou-Muratorio, K.-H. Rexer, M. Rudawska, M. Schädler, G. Taylor, K. Tuba, M. Viger, F. Villani, M. Villar</footnote>
First, these genetic changes and changes in the associated
traits may lead to shifts in the occurrence and abundance
of species already associated with the host. Second, the
new phenotypic traits of the focal species may allow new
species from the regional species pool to colonize it.
Finally, changes in the genetics of the host may induce
evolutionary responses, including speciation events, in the
associated organisms, which may feedback to evolutionary
changes in the host.

If the above scenarios hold true, we expect the relatedness
of host genotypes to co-vary with similarity among the
communities of associated species (Bangert et al. 2006;
Brändle & Brandl 2006). Within species, such patterns

have received considerable attention under the concept
of the “extended phenotype”. This concept was
introduced by Richard Dawkins (1982) to describe effects
of genes on an individual’s environment including other
organisms. Whitham et al. (2003; 2005; 2006) adopted
this concept and developed a framework for community
and ecosystem genetics, which includes a feedback where
an individual’s phenotype is dependent on the interaction
with other species.

Community assembly (Kraft et al. 2007; Emerson & Gillespie 2008) is shaped by successive filters, including regional species pool, habitat area and isolation (biogeographical filters), local environmental constraints

<caption>FIGURE 1</caption>

How host plant genes might shape assemblages of associated organisms (blue pathway on the left). Several ecological filters drive the structure of communities associated with one host plant. Among associated species co-occurring within a region and determined by evolutionary and biogeographical processes (1, Total species pool), local species assemblages depend on dispersal (2, Landscape species pool) and habitat filters (3, Habitat species pool). Dispersal filter refers to the ability of species to colonize the focal site. Habitat filters correspond to their capacity to develop and survive in a habitat given abiotic constraints. Biotic interactions with the host species contribute to the shaping of a host species pool (4, biotic filter). Finally, variation among host plant genotypes may further select different associated communities, shaping the extended phenotypes.

Genes of the focal host plant can interact with the four
filters, as illustrated by the interaction between trees
and associated insect herbivores: (1) There is evidence

that pools of insect herbivore species of different tree
families or genera are significantly different, probably
owing to a long co-evolutionary process involving
insect feeding traits and plant defence responses
(Novotny et al. 2002); (2) insect herbivores use
genetically controlled physical (e.g. shape, colour) and
chemical cues (e.g. volatile organic compounds)
provided by host plants to locate the plants; (3) trees
can be seen as ecological engineers which can modify
abiotic conditions that insects experience, e.g. wind,
moisture, or light; (4) genes control plant phenotype
and resistance traits that are deeply involved in
interactions with insect herbivores (Schoonhoven et
al. 2005); and (5) variants of host plant genes may
ultimately induce quantitative changes in traits
involved in plant-insect interactions with consequences
for insect community structure (Crutsinger et al.
2008).

Presumed reciprocal effects, through which
associated organisms feed back to the composition of
host genes, are depicted by orange colors (right side).

<figure><img src="image_66.png" /></figure>
(abiotic, biophysical filters) and biological interactions such as competition or predation (biotic filters; Fig. 1). The host genotype, interacting with the environment, may affect the structure of associated communities at several filtering steps by controlling phenotypic traits that allow associated organisms to locate, select and exploit resources of their host (Johnson & Agrawal 2005; Bailey et al. 2009) (Fig. 1). Thus, spatial variation in the composition of associated communities has a strong regional component.

Despite many reports demonstrating a correlation between genotypes of a focal species and the composition of associated communities, the fundamental ecological, genetic and evolutionary processes that generate this correlation remain poorly explored and require consideration in future studies. In this regard, three aspects deserve special attention: spatial variation, temporal variation, and gene-to-gene interactions.

First, space needs to be better integrated into study designs. As noted above, the assembly of species depends on the regional species pool, whose phylogenetic and functional structure imposes a constraint on the emerging local communities (Fig. 1). A group of genotypes of a focal species in natural or experimental populations is embedded in a landscape context that may include forest patches, arable land, urban environment or other habitats type, each of which has different species pools that might interact with the focal species. As the associated community influences the fitness of the focal species, the relative fitness of these genotypes will vary across sites, even if the abiotic conditions are similar. However, in single common garden experiments, genotypes of a focal species are exposed only to one particular species pool. Therefore, regional replicates of such experiments are necessary to estimate the stability of relationships between genotypes of the focal species and communities of associated species. Such replicates would enable us to distinguish between mainly spatial effects and those that can be attributed to the interaction between host genotypes and associated organisms. Alternatively, one might set up more complex common gardens including particular treatments, for example through fertilization or irrigation. Such an approach would allow tests of the effect of genotype x environment interactions on the assemblage of associated species for each local species pool. Furthermore, replicated common garden experiments would allow constructing reaction norms of different genotypes of the focal species. Do these genotypes respond differentially for their extended phenotypes to the changes of abiotic or biotic conditions across the testing sites? An initial step would be to identify the shape of the reaction norms (linear or quadratic) and then to estimate their variation among genotypes. Finally, the spatial context may also be dissected at the within-population level. For example, natural populations of trees usually exhibit strong spatial autocorrelation due to limited

dispersal, which increases steadily over generations. On the other hand random spatial genetic structure is observed in recently planted forests. One would therefore expect very different spatial structures of extended phenotypes among these strongly contrasting cultural regimes.

Second, community genetics should consider temporal variation in species interactions, e.g. among seasons, among years along successional sequences, and other types of temporal gradients. Traits involved in plant-herbivore interactions are known to change during plant ontogeny (Boege & Marquis 2005; Holeski et al. 2009), which is why communities of insect herbivores - and herbivory pressure - on seedlings and mature individuals may differ (Le Corff & Marquis 1999; Basset 2001). Furthermore, although associated communities may change within and between years due to fluctuations in plant phenotypes, equally they may change due to differences in weather conditions. Thus, the phenotypic traits that are important for species interactions in a particular season or year may change within and between years, and drawing conclusions from short-term experiments may be misleading. Although such traits, and the underlying genes, are genuinely involved in community interactions, their relative importance compared to other genes may vary in time and can therefore only be established in long-term experiments. Hundreds of insect generations interact with a long-living host such as a tree during its lifetime, and each generation experiences different biophysical constraints and trophic interactions with other fungi, herbivores or predators. As a consequence, even though insect populations can adapt to individual host genotypes (Mopper et al. 2000), the strength and direction of these adaptations are likely to change over time (moving targets; Ruhnke et al. 2006).

Moreover, genetic processes underlie the formation of adaptive demes and co-evolution between host and associated organisms (Fig. 1). At present, the number and type of genes involved and the associated phenotypes of interacting species are largely unknown. Recent technological advances enable researchers to sequence whole genomes and to monitor gene expression of interacting species, offering the potential to identify the candidate genes mediating the interactions between focal and associated species. Such approaches will move community genetics from studying anonymous genotype/phenotype effects to studying gene-to-organism, gene-to-gene, and ultimately to genome-to-genome interactions. While current research has focused on the few "genome-enabled" species (Ekblom & Galindo 2011), the many ongoing whole-genome projects will widen the array of study systems applying genomics data in the near future (e.g. http://www.arthropodgenomes.org/wiki/i5K, http://1000.fungalgenomes.org/home/, http://pinegenome.org/pinerefseq/).
The following sections describe how the various types of -omics may stimulate community genetics. and how they enable the genetic component of variation in community composition to be addressed at the level of variants in adaptive genes and their differential expression.

## An example of functional genomics based on complete genome sequences: ectomycorrhizal symbiosis

Ectomycorrhizae, the mutualistic symbiosis between tree roots and a cortege of soil fungal partners, are the most widespread and species-rich associations in temperate and boreal forests. Ectomycorrhizal fungi receive carbon from photosynthesis and, in turn, promote tree growth, enhance the survival of seedlings and increase the fitness of their plant partners under a wide range of environmental conditions. Despite the ecological significance of this mutualistic interaction, we have only started to explore its role for community ecology.

A breakthrough was the release of the first two full-genome drafts of mycorrhizal fungi, namely *Laccaria bicolor* (Basidiomycota) and *Tuber melanosporum*, the Périgord truffle (Ascomycota; Martin et al. 2008; Martin et al. 2010). Comparative genomics of the two mycorrhizal fungi indicated that they use different gene networks ('molecular toolkits') to establish symbiosis (Martin et al. 2010). There are vast differences between these two ectomycorrhizal genomes. *Laccaria bicolor* has a 65 Mb genome with more than 23 000 predicted proteins, which is the largest complement of genes known for any fungus, whereas *T. melanosporum* has the largest fungal genome so far with 125 Mb, but has only 7500 predicted genes, one of the smallest complement of proteins in any filamentous fungal genomes sequenced so far. Also, whereas the secretion of effector-like small secreted proteins seems to be crucial for the establishment of the symbiosis in *L. bicolor* (Plett et al. 2011), these so-called mycorrhizal-induced small secreted proteins (MiSSPs) are not present in the transcriptome of *T. melanosporum* symbiotic tissue (Martin et al. 2010). In spite of these differences, some common features and some novelties emerged from the comparison with genomes of saprophytic and pathogenic fungi. Besides the loss of plant cell-wall degrading enzymes in ectomycorrhizae, an increase in the diversity and expression of nutrient transporters and signalling pathways (e.g. tyrosine kinases) in symbiotic tissues are hallmarks of mycorrhizal genomes (Martin et al. 2008; Kosti et al. 2010; Martin et al. 2010; Plett et al. 2011). These symbiosis-related genes are good candidates for gene expression studies of multi-species interactions in the field. On the tree side, it is not known how the host tree selects its symbiotic associates. Plant-encoded small secreted proteins may be required, as shown for nitrogen-fixing symbioses (Van de Velde et al. 2010). Genomic studies will probably be the only way to elucidate the mechanisms of

interaction and to understand the effect of gene variants on this interplay. Therefore, we think that this system is an exciting model for community genetics in the -omics era.

Ectomycorrhizal fungi show a continuum of specialization to the host tree from strict specialists to generalists. Differences in the expansion of multigene families, in particular dynamic repertoires of genes encoding small secreted proteins and sugar-cleaving enzymes, might be responsible for the different host ranges of specialists, such as *T. melanosporum*, and generalists, such as *L. bicolor* (Martin et al. 2010). That is, the genome expansion observed in *L. bicolor* might be driven by selection of the symbiont to exploit diverse substrates provided by multiple potential hosts and by diverse soils. As more genomes of mycorrhizal fungi are sequenced (Martin et al. 2011), this hypothesis will become testable.

In addition to the genomics of host-symbiont interactions, studies of geographical patterns of co-evolution add to our knowledge of processes leading to reciprocal adaptation and specialization. There are only a handful of studies reporting the structure of geographic variation and patterns of co-evolution in mycorrhizal interactions, indicating that these patterns are geographically highly variable (Hoeksema 2010; Hoeksema et al. 2012). To date, mostly higher-level traits, such as intensity of mycorrhizal colonization or growth of host trees, have been studied. Several of these studies found significant genetic variation in either the host plant or the mycorrhizal fungus in its ecological effect on the other partner. For example, the relationship between the colonization intensity of the ectomycorrhizal fungus *Thelephora terrestris* and the growth of its host, Lodgepole pine (*Pinus contorta*), depends on the tree's genotype (Karst et al. 2009). In poplar, both the intensity of colonization and the amount of enzymes secreted by poplar root tips colonized by *L. bicolor* are under the genetic control of the host (Courty et al. 2011). Similar findings come from arbuscular mycorrhizal systems, where host identity has a strong effect on the fitness of different strains of *Glomus intraradices* (Ehinger et al. 2009).

An increasing body of evidence shows that subtle intraspecific differences in the genome of host plants determine the composition of interacting communities in mycorrhizal fungi (e.g. Korkama et al. 2006; Whitham et al. 2006; Sthultz et al. 2009; Karliński et al. 2010; Leski et al. 2010; Hoeksema et al. 2012). We have experimental evidence that such as intraspecific genetic variation in the host also affects the composition of interacting mycorrhizal populations (Hoeksema & Thompson 2007), but this has not yet been tested under natural conditions. To understand the links between structure and diversity of communities and ecosystem functioning, we need to know more about spatio-temporal patterns of genetic variation. There are indications that both interspecific (e.g. van der Heijden et al. 1998; Maherali & Klironomos 2007) and intraspecific (e.g. Johnson et al. 2012) diversity of
mycorrhizal fungi can regulate productivity and ecosystem functioning. We advocate studies of community and population diversity in forests and combining them with functional field studies, involving both partners of ectomycorrhizal symbioses. Numerous new techniques are emerging for gene expression studies, marker gene evaluation using comparative genomics, and enzyme activity profiling of whole ectomycorrhizal assemblages (Courty et al. 2010). The rapid development of high-throughput sequencing technologies facilitates the survey and comparison of whole microbial communities (Buée et al. 2009), although analysis, interpretation, and publication of data still needs to be optimized (Henrik Nilsson et al. 2012). Nevertheless, combined genotypic and functional studies are now feasible and may be expanded to natural and experimental gradients. Several reports indicate that soil microbe and mycorrhizal diversity differentially affect ecosystem functioning under different environmental conditions, e.g. nutrient status (van der Heijden et al. 2008). We also know that plant-associated microorganisms are an important factor influencing plant responses to climate change (Courty et al. 2010; Pickles et al. 2012). Combined genotypic and functional studies in diverse environments will help to understand current patterns and to predict changes and effects in the future.

## Associations between genes and traits: potential of next generation approaches in community genetics

An essential part of future studies in community genetics will be to identify the genes that underlie the traits of hosts that affect associated organisms. For this, sequencing of the complete genome of a host species is not sufficient. Rather, it is essential to link the presence or action of particular variants of genes or genomic regions of a host plant to the presence or abundance of associated organisms or arrays of their genes. There are basically two strategies for this, namely QTL mapping and genome-wide association studies (GWAS). We briefly outline and illustrate below the pros and cons of these two approaches for community genetics.

An example of QTL mapping of community traits of poplar is a study aimed at identifying genomic regions associated with susceptibility to insects (DeWoody et al. unpubl.data). Parents and progeny of a poplar (Populus trichocarpa × P. deltoides) F2 mapping population were assessed for various categories of leaf damage, including chewers and skeletonizers. The damage levels significantly varied among offspring genotypes. Each category was treated as a quantitative trait in a QTL mapping approach and more than ten QTLs were detected. QTLs also varied seasonally, suggesting that the insect community responds to traits and the underlying genetic variation over time. This underlines the importance of considering temporal variation in studies of community genetics, as noted above.

Another example is a study on QTLs affecting ectomycorrhizal symbiosis in a *P. deltoides* × *P. trichocarpa* F1 population (Labbé et al. 2011). Four identified QTLs were associated with candidate genes, and differential transcript levels were assessed with the help of a whole-genome microarray. The transcripts with the highest overrepresentation were, based on their gene ontology, in the repress defense mechanisms and in pathogen resistance.

Relatively few mapping populations have been produced for long-lived tree species, due to the length of time needed to maintain and study them, and the high costs associated with it. As a single cross will not contain all alleles present in a large population of an outcrossing species, not all QTLs can be detected in a single cross, and most QTL interactions will go unnoticed. Hence, several populations are necessary, and producing them would be an important investment. Next to full-sib families it may be possible to use full or partial diallel designs with multiple parents, so that more alleles are included and many more allele combinations can be studied, similar to MAGIC populations (Kover et al. 2009) but without the need for selfing to multiply and maintain the population.

In the meantime, an elegant alternative for forest trees is to use existing progeny trials. Many of these have been established and often replicated at different locations, and phenotypic data are usually available for extensive periods of time. Many trials consist of half-sib families, in which the alleles from the mother segregate in the progeny. If only a limited number of fathers were involved, genotyping may even allow them to be split into a few interconnected full-sib families. Common garden experiments often include a sample of the diversity of an area. When these experiments are replicated at multiple sites, it may be possible to perform genome-wide association mapping with the advantage of multi-site / multi-year data.
An issue for community genetics, as mentioned above, is that the local species pool may be different between the locations of the trials. This can be tackled efficiently by replicating the populations and planting them in different locations. Replicated populations will also spread the risk of losing individual members of the populations.

After finding a QTL region based on the presence of an associated organism or, for example, damage caused by an insect species, the underlying mechanism can be unravelled, in this case by measuring the secondary compound composition of all progeny trees and locating such traits on the genetic map. Co-localization of a compound with a QTL would suggest that it was responsible for the effect on the insects and that a structural or regulatory gene involved in its synthesis is located in that genomic region. In some species, this can be tested by mutant analysis, but it is not practical with trees. Alternatively, one could analyse the naturally occurring genetic variation in a large set of unrelated trees
<caption>TABLE 1</caption>

From genes of focal species to traits of the extended phenotype – and back: questions and experimental considerations, related to (a) spatio-temporal variation, (b) the application of -omics approaches, and (c) reciprocal effects to stimulate future studies in community genetics

<table><thead><tr><th>Theme</th><th>Questions</th><th>Experimental considerations</th></tr></thead><tbody><tr><td rowspan="5">(a) Spatio-temporal variation</td><td>To what degree do regional species pools determine the composition of organisms associated to particular genotypes?</td><td>Assess naturally occurring spatial replicates of particular genotypes, e.g. agricultural, horticultural or silvicultural clones, and perform regionally replicated experiments using the same (set of) genotypes exposed to various regional species pools of potentially associated organisms.</td></tr><tr><td>What is the relevance of phylogeographic structure in host species for the composition of associated communities?</td><td>Consider genetic structure and evolutionary lineages of the focal species.</td></tr><tr><td>How do relationships between genotypes and associated organisms vary among seasons or among life stages?</td><td>Perform temporally replicated experiments or monitor natural communities across >1 year; establish long-term experiments with host plants from seedlings to mature adults.</td></tr><tr><td>How does landscape configuration, e.g. differences in the relative abundance of, or connectivity among, particular habitat types, affect regional species pools and, thus, the communities of associated organisms in a focal species?</td><td>Include landscape characteristics when setting up experimental plots or assessing natural communities.</td></tr><tr><td>To what degree does phenotypic plasticity shape extended phenotypes?</td><td>Set up common garden experiments along ecological gradients including reciprocal transplants to test for genotype-by-environment interactions and reaction norms.</td></tr><tr><td rowspan="4">(b) -omics approaches</td><td>Which QTL relate to particular groups of associated organisms?</td><td>Establish various fullsib families or diallel crosses to include a wide range of allele variants.</td></tr><tr><td>What (classes of) compounds differ among host genotypes that are differentially affected by groups of associated organisms?</td><td>Genome/transcriptome sequencing of pools of host plants differing in their associated communities.</td></tr><tr><td>Do traits affecting community composition of associated species rely on single or multiple genes, and how large is their allelic variation within host populations?</td><td>Identify genes directly involved in the interaction, e.g. through QTL mapping, and quantify the degree of polymorphism using high-throughput, reduced-representation sequencing.</td></tr><tr><td>Does one gene of a focal species influence a single, a group of, or all associated species?</td><td>Use feeding (herbivores) or inoculation (ectomycorrhizae) experiments and perform co-expression profiling and subsequent protein annotation.</td></tr><tr><td rowspan="3">(c) Reciprocal effects</td><td>How many such genes exist, given that a focal species may interact with hundreds of associated species?</td><td>Perform gene expression studies of focal species that are experimentally associated with different single species or groups of species of associated organisms.</td></tr><tr><td>How do different groups of associated species induce changes in the phenotypic traits (and the underlying allele frequencies) of the host?</td><td>Expose the same (set of) hosts to different (sets of) associated species and test for changes in traits and allele frequencies over time.</td></tr><tr><td>What genes in host and associated species determine whether they interact as generalists or specialists?</td><td>Combine comparative genomics and expression profiling among generalists and specialists in both hosts and associated species.</td></tr></tbody></table>
with different combinations of compounds and conduct association tests (i.e. GWAS).

Genome-wide association studies assume that, in the absence of population substructure, markers that are physically linked to a gene associated with a phenotype of a trait can be distinguished from markers that are not linked, as the latter are assumed to occur randomly in individuals of the population regardless of the phenotype (Nordborg & Weigel 2008). There is no need to construct a mapping population as in QTL detection, but a reference genome or a dense genetic map in combination with sufficient linkage disequilibrium (LD) are required (Kim et al. 2007). LD appears to be limited in tree species (Ingvarsson 2005; Heuertz et al. 2006; Pyhäjärvi et al. 2007), which implies that high-density genetic marker arrays are needed for applying association mapping and that many more individuals need to be studied. For instance, Fournier-Level et al. (2009) tested target candidate genes and identified the functional variation responsible for the observed variation in anthocyanin variation in grape by association analysis. The very low LD often encountered in natural tree populations (Neale & Savolainen 2004) will assist in finding many of the possible combinations of compounds, thus increasing the power of the association study. A new approach, becoming feasible because of high-throughput sequencing technology, is to pool and sequence DNA from multiple individuals within a population with clearly distinct phenotypes or habitat conditions (Turner et al. 2010), and to identify those markers across the genome that display a large difference in allelic frequency between the pooled groups (Holderegger et al. 2008). The advantage of this 'population resequencing' approach, which vaguely resembles bulked segregant analysis (BSA), is that no mapping population or extensive LD is necessary; the drawback is that an annotated genome is still needed for reference. Since annotated genome sequences are increasingly becoming available, this will be less of a problem in the future. The approach can be readily extended to polygenic traits (Heard et al. 2010). A potential application to community genetics in trees would be to pool the DNA from trees that host a particular insect with DNA from those that do not, and compare the sequenced genomes of the two groups.

Next generation methods now enable genotyping-by-sequencing (Baird et al. 2008). In the context of segregating populations, restriction-site associated DNA (RAD) markers or transcriptome sequencing permit direct mapping-by-sequencing, thus skipping marker development altogether (Hartwig et al. 2012; Zhu et al. 2012). In QTL mapping this solves the problem of generating dense maps, so that the limiting factor for high resolution is the number of recombinations or the size of the segregating population. As forest trees have very small LD, the ability to generate high volumes of genomic data is a very promising development for GWAS.

Gene expression profiling, a complementary approach to association genomics as a strategy for functional genomics, is also being revolutionized by developments in next generation technologies. Gene expression profiling has been applied to study stress response in trees, for example following insect attack where transcript analyses by cDNA microarray profiles have been combined with 2-D protein and protein spectrometric analyses (Lippert et al. 2007). In this pioneering work on pines and pine weevils, the authors demonstrated that transcripts and their proteins were complimentary. Next generation sequencing of tagged cDNA ends now enables researchers to quantify the number of transcripts from different subsets of individuals (Xu et al. 2009). Given the availability of gene annotations, the transcripts will be associated with gene models and their regulators using publicly available databases. We expect that co-expression profiling will become feasible for populations as well as for individual ontogenetic stages of interacting species. Such an approach may also be scaled up from two-species interactions to multiple-species interactions, i.e. a true 'community transcriptome' approach.

Proteomic approaches allow for an efficient and simultaneous detection of the proteins in a sample. The proteome composition to some extent integrates fluctuations in expression over a period of time, thus potentially being robust with regard to sampling time in the field. The identification of peptides relies on either a large, high-quality RNA-seq dataset, a complete set of alleles from a multigene family, or the genome sequence. An example is the use of peptide identification (Q-TOF LC-MS<sup>E</sup>) for fast screening of Bet v 1 isoforms in pollen of various birch species, as it was possible to determine both presence and relative abundances of individual isoforms (Schenk et al. 2009). For this, the mass spectra obtained from the pollen were compared with a set of predicted peaks based on a complete set of isoforms obtained by sequencing the genes. In species for which the genome sequence or a large amount of transcriptome data is available, this prediction becomes a relatively simple bioinformatics exercise.

Other omics techniques, such as metabolomics, may be employed in similar experimental schemes. Recent advances have increased the sensitivity and throughput of metabolomics and proteomics assays ('next-gen biochem'). Now, one can directly map QTL controlling the metabolic profile of all offspring of a cross. For instance, untargeted GC-TOF-MS metabolite profiling allowed mapping of 100 mQTLs (Carreno-Quintero et al. 2012). The main drawbacks of metabolomics are the higher costs and the problem of interfering factors due to the different growing conditions of the trees included in the association analysis. Moreover the samples cannot be all taken at the same time. On the other hand, the compounds measured are also the ones that affect the interaction with associated insect species. So if genetic variation in multiple
genes affects the content of one important compound, the association of the compound with presence or absence of one or more insect species will be stronger than that of each of the underlying genes, and the association will also be more informative on the mechanism of the interaction. Even GWAS could be done in this way. In our example using a pool of trees including those that host a particular insect and those that do not, a comparison of compounds may be more straightforward than comparing DNA markers. In particular, if the insect is not always present on the same trees across years, the compounds present in each tree in each year could reveal a strong correlation, whereas the genes that enable the tree to produce the compounds would not.

If, as indicated above, a compound affects the presence of insect species, then one would expect, reciprocally, the presence of catabolites of the compound in insect species that tolerate the compound, when these insects are sampled on the trees that produce it. This can be used to experimentally validate the statistical associations between compounds in the tree and the presence of insect species or guilds, and for a starting point for understanding the mechanisms behind the interactions between trees and insects.

## Perspectives

A suite of -omics approaches is available to pave the way for studying entire communities. Accordingly, we need to refine hypotheses and develop suitable study designs and statistical tools (Augustin et al. 2010; Ovaskainen et al. 2010), which will improve implementation once reduced costs make these tools applicable to large-scale sampling of community-level interactions (Table 1).

As outlined above, we see two main directions that should be followed in community genetics to substantiate inference on the interplay of genes, organisms, communities, and their respective environments. First, joint descriptive and experimental studies should include spatial and temporal gradients to account for environmental variation in these dimensions (Thompson 2005; Crutsinger et al. 2009; Tack et al. 2010). Second, researchers in community genetics should make better use of the exponentially increasing genomic information becoming available, which will requires solid expertise in bioinformatics. If this is achieved, gene-to-gene interactions can be explored in individual-based associations and at the level of entire communities and shift community genetics towards becoming community genomics.

Moreover, community genetics goes beyond the effects of genotypes in one species on the community of associated organisms. We also need to consider the reciprocal effects of how associated communities shape the genotypic composition of their hosts and of how the genotypes of associated species affect host communities (Fig. 1). There

are virtually no studies available on this aspect of community interactions, which leaves a wide-open field of empirical research for the future. Exploring reciprocal interactions might help to extrapolate population genomics and quantitative genomics of focal species. We will then need to adopt a community-based understanding of selection and drift as well as to include G x G x E interactions into reaction norm calculations. However, elaborating on this subject goes beyond the scope of the present article.

In conclusion, we believe that the amalgamation of traditional population genetics, quantitative genetics and ecology, fostered by the advent of new genomic technologies, will revolutionise our perception of community and ecosystem processes and push community genetics into a new era.

## Acknowledgments

Our sincere thanks go to Rolf Holderegger and the anonymous reviewers for helpful comments on earlier versions of the manuscript. Karen A. Brune and Stephen Cavers provided linguistic editing. This work was financed through the EC-supported Network of Excellence Evoltree (GOCE-016322).

## REFERENCES

Agrawal AA (2003) Community genetics: new insights into community ecology by integrating population genetics. *Ecology*, **84**, 543-544.

Augustin S, Denux O, Castagneyrol B et al. (2010) Insect herbivory response to Populus nigra genetic diversity. Fifth International Symposium, Poplars and willows: from research models to multipurpose trees for a bio-based society. Orvieto, Sept. 20-25, 2010, p. 207.

Bailey JK, Hendry AP, Kinnison MT et al. (2009) From genes to ecosystems: an emerging synthesis of evolutionary dynamics. *New Phytologist*, **184**, 746-749.

Baird NA, Etter PD, Atwood TS et al. (2008) Rapid SNP discovery and genetic mapping using sequenced RAD markers. *PLoS One*, **3**, e3376.

Bangert RK, Turek RJ, Rehill B, Wimp GM, Bailey JK, Whitham TG (2006) A genetic similarity rule determines arthropod community structure. *Molecular Ecology*, **15**, 1379-1391.

Basset Y (2001) Communities of insect herbivores foraging on saplings versus mature trees of *Pourouma bicolor* (Cecropiaceae) in Panama. *Oecologia*, **129**, 253-260.

Boege K, Marquis RJ (2005) Facing herbivory as you grow up: the ontogeny of resistance in plants. *Trends in Ecology & Evolution*, **20**, 441-448.

Brändle M, Brandl R (2006) Is the composition of phytophagous insects and parasitic fungi among trees predictable? *Oikos*, **113**, 296-304.

Buée M, Reich M, Murat C et al. (2009) 454 Pyrosequencing analyses of forest soils reveal an unexpectedly high fungal diversity. *New Phytologist*, **184**, 449-456.
Carreno-Quintero N, Acharjee A, Maliepaard C et al. (2012) Untargeted metabolic quantitative trait loci analyses reveal a relationship between primary metabolism and potato tuber quality. *Plant Physiology*, **158**, 1306-1318.

Courty P-E, Buée M, Diedhiou AG et al. (2010) The role of ectomycorrhizal communities in forest ecosystem processes: New perspectives and emerging concepts. *Soil Biology & Biochemistry*, **42**, 679-698.

Courty PE, Labbé J, Kohler A et al. (2011) Effect of poplar genotypes on mycorrhizal infection and secreted enzyme activities in mycorrhizal and non-mycorrhizal roots. *Journal of Experimental Botany*, **62**, 249-260.

Crutsinger GM, Cadotte MW, Sanders NJ (2009) Plant genetics shapes inquiline community structure across spatial scales. *Ecology Letters*, **12**, 285-292.

Crutsinger GM, Reynolds WN, Classen AT, Sanders NJ (2008) Disparate effects of plant genotypic diversity on foliage and litter arthropod communities. *Oecologia*, **158**, 65-75.

Dawkins R (1982) The Extended Phenotype: *The Gene as the Unit of Selection*. Freeman, Oxford.

Ehinger M, Koch AM, Sanders IR (2009) Changes in arbuscular mycorrhizal fungal phenotypes and genotypes in response to plant species identity and phosphorus concentration. *New Phytologist*, **184**, 412-423.

Ekblom R, Galindo J (2011) Applications of next generation sequencing in molecular ecology of non-model organisms. *Heredity*, **107**, 1-15.

Emerson BC, Gillespie RG (2008) Phylogenetic analysis of community assembly and structure over space and time. *Trends in Ecology & Evolution*, **23**, 619-630.

Fournier-Level A, Cunff LL, Gomez C et al. (2009) Quantitative genetic bases of anthocyanin variation in grape (*Vitis vinifera* L. ssp. *sativa*) berry: A quantitative trait locus to quantitative trait nucleotide integrated study. *Genetics*, **183**, 1127-1139.

Hartwig B, Velikkakam James G, Konrad KK, Schneeberger K, Turck F (2012) Fast isogenic mapping-by-sequencing of ethyl methanesulfonate-induced mutant bulks. *Plant Physiology*, **160**, 591-600.

Heard E, Tishkoff S, Todd JA et al. (2010) Ten years of genetics and genomics: what have we achieved and where are we heading? *Nature Reviews Genetics*, **11**, 723-733.

Henrik Nilsson R, Tedersoo L, Lindahl BD et al. (2012) Towards standardization of the description and publication of next-generation sequencing datasets of fungal communities. *New Phytologist*, **191**, 314-318.

Heuertz M, de Paoli E, Källmann T et al. (2006) Multilocus patterns of nucleotide diversity, linkage disequilibrium and demographic history of Norway spruce [*Picea abies* (L.) Karst]. *Genetics*, **174**, 2095-2105.

Hoeksema JD (2010) Ongoing coevolution in mycorrhizal interactions. *New Phytologist*, **187**, 286-300.

Hoeksema JD, Hernandez JV, Rogers DL, Mendoza LL, Thompson JN (2012) Geographic divergence in a species-rich symbiosis: Interactions between Monterey pines and ectomycorrhizal fungi. *Ecology*, **93**, 2274-2285.

Hoeksema JD, Thompson JN (2007) Geographic structure in a widespread plant-mycorrhizal interaction: pines and false truffles. *Journal of Evolutionary Biology*, **20**, 1148-1163.

Holderegger R, Herrmann D, Poncet B et al. (2008) Land ahead: using genome scans to identify molecular markers of adaptive relevance. *Plant Ecology and Diversity*, **1**, 273-283.

Holeski LM, Kearsley MJC, Whitham TG (2009) Separating ontogenetic and environmental determination of resistance to herbivory in cottonwood. *Ecology*, **90**, 2969-2973.

Ingvarsson PK (2005) Nucleotide polymorphism and linkage disequilibrium within and among natural populations of European aspen (*Populus tremula* L., Salicaceae). *Genetics*, **169**, 945-953.

Johnson D, Martin F, Cairney J, Anderson IC (2012) The importance of individuals: intraspecific diversity of mycorrhizal plants and fungi in ecosystems. *New Phytologist*, **194**, 614-628.

Johnson MTJ, Agrawal AA (2005) Plant genotype and environment interaction to shape a diverse arthropod community on evening primrose (*Oenothera biennis*). *Ecology*, **86**, 874-885.

Johnson MTJ, Stinchcombe JR (2007) An emerging synthesis between community ecology and evolutionary biology. *Trends in Ecology & Evolution*, **22**, 250-257.

Karliński L, Rudawska M, Kieliszewska-Rokicka B, Leski T (2010) Relationship between genotype and soil environment during colonization of poplar roots by mycorrhizal and endophytic fungi. *Mycorrhiza*, **20**, 315-324.

Karst J, Jones M, Turkington R (2009) Ectomycorrhizal colonization and intraspecific variation in growth responses of lodgepole pine. *Plant Ecology*, **200**, 161-165.

Kim S, Plagnol V, Hu TT et al. (2007) Recombination and linkage disequilibrium in *Arabidopsis thaliana*. *Nature Genetics*, **39**, 1151-1155.

Korkama T, Pakkanen A, Pennanen T (2006) Ectomycorrhizal community structure varies among Norway spruce (*Picea abies*) clones. *New Phytologist*, **171**, 815-824.

Kosti I, Mandel-Gutfreund Y, Glaser F, Horwitz B (2010) Comparative analysis of fungal protein kinases and associated domains. *BMC Genomics*, **11**, 133.

Kover PX, Valdar W, Trakalo J et al. (2009) A multiparent advanced generation nter-cross to fine-map quantitative traits in *Arabidopsis thaliana*. *PLoS Genetics*, **5**, e1000551.

Kraft NJB, Cornwell WK, Webb CO, Ackerly DD (2007) Trait evolution, community assembly, and the phylogenetic structure of ecological communities. *American Naturalist*, **170**, 271-283.

Labbé J, Jorge V, Kohler A et al. (2011) Identification of quantitative trait loci affecting ectomycorrhizal symbiosis in an interspecific F1 poplar cross and differential expression of genes in ectomycorrhizas of the two parents: *Populus deltoides* and *Populus trichocarpa*. *Tree Genetics & Genomes*, **7**, 617-627.

Le Corff J, Marquis RJ (1999) Differences between understorey and canopy in herbivore community composition and leaf quality for two oak species in Missouri. *Ecological Entomology*, **24**, 46-58.

Leski T, Aučina A, Skridaila A, Pietras M, Riepšas E, Rudawska M (2010) Ectomycorrhizal community structure of different genotypes of Scots pine under forest nursery conditions. *Mycorrhiza*, **20**, 473-481.
Lippert D, Chowrira S, Ralph SG et al. (2007) Conifer defense against insects: Proteome analysis of Sitka spruce (*Picea sitchensis*) bark induced by mechanical wounding or feeding by white pine weevils (*Pissodes strobi*). *Proteomics*, **7**, 248-270.

Maherali H, Klironomos JN (2007) Influence of phylogeny on fungal community assembly and ecosystem functioning. *Science*, **316**, 1746-1748.

Martin F, Aerts A, Ahren D et al. (2008) The genome of *Laccaria bicolor* provides insights into mycorrhizal symbiosis. *Nature*, **452**, 88-92.

Martin F, Cullen D, Hibbett D et al. (2011) Sequencing the fungal tree of life. *New Phytologist*, **190**, 818-821.

Martin F, Kohler A, Murat C et al. (2010) Perigord black truffle genome uncovers evolutionary origins and mechanisms of symbiosis. *Nature*, **464**, 1033-1038.

Mopper S, Stiling P, Landau K, Simberloff D, Van Zandt P (2000) Spatiotemporal variation in leafminer population structure and adaptation to individual oak trees. *Ecology*, **81**, 1577-1587.

Nadeau NJ, Jiggins CD (2010) A golden age for evolutionary genetics? Genomic studies of adaptation in natural populations. *Trends in Genetics*, **26**, 484-492.

Neale DB, Savolainen O (2004) Association genetics of complex traits in conifers. *Trends in Plant Science*, **9**, 325-330.

Nordborg M, Weigel D (2008) Next-generation genetics in plants. *Nature*, **456**, 720-723.

Novotny V, Basset Y, Miller SE et al. (2002) Low host specificity of herbivorous insects in a tropical forest. *Nature*, **416**, 841-844.

Ovaskainen O, Nokso-Koivista J, Hottola J et al. (2010) Identifying wood-inhabiting fungi with 454 sequencing - what is the probability that BLAST gives the correct species? *Fungal Ecology*, **3**, 274-283.

Pickles BJ, Egger KN, Massicotte HB, Green DS (2012) Ectomycorrhizas and climate change. *Fungal Ecology*, **5**, 73-84.

Plett J, Kemppainen M, Kale S et al. (2011) A secreted effector protein of *Laccaria bicolor* is required for symbiosis development. *Current Biology*, **21**, 1197-1203.

Pyhäjärvi T, Rosario García-Gil M, Knürr T, Mikkonen M, Wachowiak W, Savolainen O (2007) Demographic history has influenced nucleotide diversity in European *Pinus sylvestris* populations. *Genetics*, **177**, 1713-1724.

Rokas A, Abbot P (2009) Harnessing genomics for evolutionary insights. *Trends in Ecology and Evolution*, **24**, 192-200.

Rowntree JK, Shuker DM, Preziosi RF (2011) Forward from the crossroads of ecology and evolution. *Philosophical Transactions of the Royal Society B: Biological Sciences*, **366**, 1322-1328.

Ruhnke H, Schädler M, Matthies D, Klotz S, Brandl R (2006) Are sawflies adapted to individual host trees? A test of the adaptive deme formation hypothesis. *Evolutionary Ecology Research*, **8**, 1039-1048.

Schenk MF, Cordewener JHG, America AHP, van't Westende WPC, Smulders MJM, Gilissen LJWJ (2009) Characterization of PR-10 genes from eight *Betula* species and detection of Bet v 1 isoforms in birch pollen. *BMC Plant Biology*, **29**, 24.

Schoonhoven LM, Loon JJAv, Dicke M (2005) *Insect-Plant Biology*, 2nd edn. Oxford University Press, New York.

Sthultz CM, Whitham TG, Kennedy K, Deckert R, Gehring CA (2009) Genetically based susceptibility to herbivory influences the ectomycorrhizal fungal communities of a foundation tree species. *New Phytologist*, **184**, 657-667.

Tack AJM, Ovaskainen O, Pulkkinen P, Roslin T (2010) Spatial location dominates over host plant genotype in structuring an herbivore community. *Ecology*, **91**, 2660-2672.

Thompson JN (2005) *The geographic mosaic of coevolution*. University of Chicago Press, Chicago.

Turner TL, Bourne EC, Von Wettberg EJ, Hu TT, Nuzhedin SV (2010) Population resequencing reveals local adaptation of *Arabidopsis lyrata* to serpentine soils. *Nature Genetics*, **42**, 260-263.

Van de Velde W, Zehirov G, Szatmari A et al. (2010) Plant peptides govern terminal differentiation of bacteria in symbiosis. *Science*, **327**, 1122-1126.

van der Heijden MGA, Bardgett RD, Van Straalen NM (2008) The unseen majority: soil microbes as drivers of plant diversity and productivity in terrestrial ecosystems. *Ecology Letters*, **11**, 296-310.

van der Heijden MGA, Klironomos JN, Ursic M et al. (1998) Mycorrhizal fungal diversity determines plant biodiversity, ecosystem variability and productivity. *Nature*, **396**, 69-72.

Whitham TG, Bailey JK, Schweitzer JA et al. (2006) A framework for community and ecosystem genetics: from genes to ecosystems. *Nature Reviews Genetics*, **7**, 510-523.

Whitham TG, Lonsdorf EV, Schweitzer JA et al. (2005) "All effects of a gene on the world": Extended phenotypes, feedbacks, and multi-level selection. *Ecoscience*, **12**, 5-7.

Whitham TG, Young WP, Martinsen GD et al. (2003) Community and ecosystem genetics: a consequence of the extended phenotype. *Ecology*, **84**, 559-573.

Wymore AS, Keeley ATH, Yturralde KM, Schroer ML, Propper CR, Whitham TG (2011) Genes to ecosystems: exploring the frontiers of ecology with one of the smallest biological units. *New Phytologist*, **191**, 19-36.

Xu WJ, Wang ZX, Qiao ZD (2009) Modified PCR methods for 3' end amplification from serial analysis of gene expression (SAGE) tags. *FEBS Journal*, **276**, 2657-2668.

Zhu Y, Mang H-g, Sun Q, Qian J, Hipps A, Hua J (2012) Gene discovery using mutagen-induced polymorphisms and deep sequencing: application to plant disease resistance *Genetics*, **192**, 139-146.

Author contribution

All authors have conceptually contributed to the content of the present article through their active involvement and the many discussions held during the EvolTree Jointly Executed Research Activity 3 (JERA 3) on community ecology, which was managed by BZ. FG designed and organized the article preparation, FG, RB, BC, HJ, MP, MJMS, and BZ drafted the various sections, and all main authors participated in the writing of the final text.
Antoine Kremer<fnref n="1" />,<fnref n="2" />*, Ophélie Ronce<fnref n="3" />*, Juan J. Robledo-Arnuncio<fnref n="4" />, Frédéric Guillaume<fnref n="5" />, Gil Bohrer<fnref n="6" />, Ran Nathan<fnref n="7" />, Jon R. Bridle<fnref n="8" />, Richard Gomulkiewicz<fnref n="9" />, Etienne K. Klein<fnref n="10" />, Kermit Ritland<fnref n="11" />, Anna Kuparinen<fnref n="12" />, Sophie Gerber<fnref n="13" />,<fnref n="2" />, Silvio Schueler<fnref n="14" />

**Corresponding Author:**

Antoine Kremer, INRA

UMR1202 Biodiversité Gènes et Communautés 69 route d'Arcachon, Cestas, F-33610, France

Email: antoine.kremer@pierroton.inra.fr

Forest trees are keystone species, dominating many terrestrial ecosystems in regions where the most pronounced climate changes will occur. Climate change generates new adaptive challenges for trees. Their life history traits could either constrain or accelerate their adaptation. On one hand, long generation time can slow down evolutionary responses. On the other hand, long distance (LD) gene flow could compensate for their long generation time, facilitating evolutionary change in a shifting climate. We critically examine the latter hypothesis, by reviewing data and theory about the extent of gene flow in trees and its evolutionary consequences. Abundant evidence of LD effective dispersal indicates that genes may move within one generation over larger scales than the predicted shifts of tree habitat. Gene flow can have antagonistic effects on adaptation and persistence in the specific temporal and spatial frame of predicted climate change. Both theory and empirical data however suggest that the positive effects of LD gene flow in forest trees may dominate in many instances. The balance between the different effects of gene flow may however differ between the leading edge, the core and the rear of the distributions. Finally, we suggest future experimental and theoretical research areas for a better integration of trees dispersal biology and evolutionary quantitative genetics.

**Key words:** selection, adaptation, gene flow, climate change, forest trees.

Published in *Ecology Letters* (2012) 15: 378-392, with doi:10.1111/j.1461-0248.2012.01746.x

Introduction

While evidence of climate change and its impact on
world's biota is steadily increasing, so are our concerns
about the biological or human mediated capacities of
species and populations to cope with these changes
(Solomon et al. 2007). Such concerns are presumably
more acute for sedentary and long lived organisms such
as trees, which are less likely to track favourable conditions
fast enough by migration. Furthermore, trees constitute a
large ecologically and economically important functional
group of woody plants that dominate many terrestrial
ecosystems in regions where the most pronounced climate
changes are projected to occur. The near-surface
temperature is expected to shift northwards in mean rates
of 110-430 m yr⁻<fnref n="1" /> during the 21st century for Mediterranean,
temperate and boreal forests, the major forest biomes in
mid and high latitudes (Loarie et al. 2009). Local estimates

of this shift, in a scale of 1 km2, vary by three orders of
magnitude (0.01 to 10 km yr⁻<fnref n="1" />). Niche modelling under
various climatic and greenhouse emission predictions
suggests that *bioclimatic envelopes* (glossary) for forest
trees will shift northwards in North America (Iverson et al.
2008) and north-eastwards in Europe (Thuillier 2003). The
estimated shift distance varies from 300 to 800 km within
one century depending on the climate or greenhouse gas
emission scenario (Mc Kenney et al. 2007), with
considerable variation across both models and species.
The large *bioclimatic envelope* of many forest trees hides
a collection of highly differentiated populations and
genotypes with contrasting adaptation to local climate
(Box 1). Shifts in *bioclimatic envelopes* are therefore likely
to generate not only potential extinction and
recolonization, but also large reorganization of genetic
diversity within the species range if divergent or locally-
adapted populations respond in different ways.

<footnote n="1">¹ INRA, UMR1202 Biodiversité Gènes et Communautés, Cestas, F-33610, France;
email: antoine.kremer@pierroton.inra.fr</footnote>

<footnote n="2">² Université de Bordeaux, UMR1202 Biodiversité Gènes et Communautés, Talence, F-33410, France;</footnote>

<footnote n="3">³ Université Montpellier 2 - CNRS, Institut des Sciences de l'Evolution, F-34095 Montpellier, Cedex 05, France; email: ophelie.ronce@univ-montp2.fr</footnote>

<footnote n="4">⁴ Department of Forest Ecology and Genetics, Forest Research Centre (CIFOR), INIA, 28040 Madrid, Spain; email: jjrobledo@gmail.com</footnote>

<footnote n="5">⁵ ETH, Department of Environmental Sciences, Universitätstrasse 16, 8092 Zürich, Switzerland; email: frederic.guillaume@env.ethz.ch</footnote>

<footnote n="6">⁶ Department of Civil, Environmental & Geodetic Engineering, Ohio State University, Columbus, OH 43210, USA; email: bohrer.17@osu.edu</footnote>

<footnote n="7">⁷ Movement Ecology Laboratory, Department of Ecology, Evolution and Behavior, The Alexander Silberman Institute of Life Sciences, The Hebrew University of Jerusalem, Edmond J. Safra Campus, Jerusalem 91904, Israel;</footnote>

<footnote>email: rnathan@cc.huji.ac.il</footnote>

<footnote n="8">⁸ School of Biological Sciences, University of Bristol, Bristol, BS8 IUG, UK; email: Jon.Bridle@bristol.ac.uk</footnote>

<footnote n="9">⁹ School of Biological Sciences and Department of Mathematics, Washington State University, Pullman, Washington 99164, USA; email: gomulki@wsu.edu</footnote>

<footnote n="10">¹⁰ INRA, UR Biostatistiques & Processus Spatiaux 546, F-84914 Avignon, France;
email: Etienne.klein@avignon.inra.fr</footnote>

<footnote n="11">¹¹ Department of Forest Sciences, University of British Columbia, Vancouver, BC V6T 1Z4, Canada; email: kermit.ritland@ubc.ca</footnote>

<footnote n="12">¹² Ecological Genetics Research Unit, Department of Biosciences, University of Helsinki, Helsinki FI-00014, Finland; e-mail: anna.kuparinen@helsinki.fi</footnote>

<footnote n="13">¹³ INRA, UMR1202 Biodiversité Gènes et Communautés, Cestas, F-33610, France; email: sophie.gerber@pierroton.inra.fr</footnote>

<footnote n="14">¹⁴ Federal Research and Training Centre for Forests, Natural Hazards and Landscape, Seckendorf-Gudent-Weg 8, 1131 Wien, Austria; email: silvio.schueler@bfw.gv.at</footnote>

<footnote>+ Co-first authors</footnote>
With climate change, environments change continuously and the optimal sets of adaptations maximizing fitness under local conditions may shift accordingly. The evolutionary responses of populations can then be pictured as a race where populations are tracking the moving optima both in time (Bürger & Krall, 2004) and in space (Pease et al. 1989; Polechova et al. 2009). Migration and adaptation are often perceived as alternative responses to these challenges (Aitken et al. 2008) because evolution allows populations to adapt to novel conditions without migrating, whereas migration lets populations track favourable conditions without evolving. Range shifts and adaptation can also occur simultaneously

(e.g. Cwynar & MacDonald 1987). Seed dispersal allows colonization of new favourable habitat. However, both seeds and pollen dispersal in trees affect the spread of genetic variation within the range. Such gene flow affects adaptation by shaping the distribution of genetic variation both within and among populations (Lopez et al. 2008). Trees are characterized by their particular life history, combining long generation time (allowing divergent strategies at different life stages) and the capacity for long dispersal distances through pollen and seeds. Given the anticipated intensity and directionality of climatic change, do trees have the adaptive capacity to respond and how will gene flow affect that response? Valuable insight into

<caption>FIGURE 1</caption>

Predicted shifts of *bioclimatic envelopes of sessile oak (Quercus petraea)* in Europe
(according to Thuiller 2003)

<figure><img src="image_67.png" /></figure>

Predicted *bioclimatic envelopes* of sessile oak in 2080, assuming that correlations between present distribution (Panel A, light grey area) and climatic data are maintained. Climate of black areas would not be suited any more to sessile oak in 2080, while climate of dark grey areas would become favourable. Overall shifts of several hundred kilometres are foreseen

which provide some hints on the scale of gene dispersal needed to track climate change. Predictions were made according to different IPCC models of greenhouse gas emissions (GG) and climatic changes (CC) (Solomon et al. 2007). Panel B: GG is A2 and CC is CSIRO2; Panel C: GG is A2 and CC is HadCM3; Panel D: GG is A1F1 and CC is HadCM3.
these issues is provided by the study of evolutionary
changes in trees during the climate change that occurred
following the last glaciations (Petit et al. 2008). Because
of their economic importance, local adaptation in forest
trees has very early on been the subject of intensive study
at very large spatial scales. Decades-old common garden
experiments of forest trees in the Northern hemisphere
suggest that an interaction between divergent selection
across contrasted environments and large pollen flow
maintained enough diversity in local forest trees
populations to support adaptability to past changing
environments (Kremer et al. 2010). Whether interaction
between gene flow and selection will be as efficient in the

future remains unclear, as the predicted rates of
environmental changes might exceed historical ones.

We here critically examine the hypothesis that long
distance (LD) gene flow could compensate for the long
generation time of trees, facilitating evolutionary change
in a shifting climate, by reviewing both data and theory,
about the extent of gene flow in forest trees and its
evolutionary consequences. We first review the recent
literature on long distance pollen and seed dispersal in
trees and show that it can match the predicted climate
change velocity within one generation. We then review the
theoretical predictions and experimental evidence for the

BOX 1

Provenance tests and norms of reaction

Provenance tests are common garden experiments that
gather usually very large number of populations planted
by forest geneticists in multiple replicates over decades
(Morgenstern 1996). They provide crucial information on
the level of genetic variation within and between
populations for fitness related traits, resulting from a
balance of divergent selection across populations, gene
flow, and random genetic drift. Extensive surveys of
genetic diversity and variation have been conducted in
these experiments and indicate that (1) extant populations
harbour large levels of genetic variation (Hamrick et al.
1992) continuously replenished by extensive gene flow (2)
adaptive traits exhibit high levels of population
differentiation, despite gene flow, as a result of strong
divergent selection (Savolainen et al. 2007) (3) clinal
patterns of population differentiation along climatic or
geographical variables are congruent across species
suggesting systematic adaptive responses to directional
selection, particularly for phenological traits; and (4) the
extant distribution of between versus within population
differentiation for fitness related traits has developed
rapidly following post glacial recolonisation and is not the
legacy of ancient population structure (Kremer et al. 2010).
Furthermore, when replicated provenance tests were
established, reaction norms of populations can be
constructed that visualise their response across a wide
range of environmental conditions (Rehfeldt et al. 1999;
Rehfeldt et al. 2002). Reaction norms of fitness-related

traits follow generally quadratic functions. Panel A
illustrates the reactions norms of two populations of *Pinus
contorta* (Rehfeldt *et al.* 1999) for height at age 20. Δ
accounts for the difference between the climate of the site
where the population stems from (dotted line) and the
optimal climate corresponding to the site where the
population exhibits the highest value for height (bold line).
Such *reaction norms* show that the climatic tolerances of
individual *provenances* are narrower than the whole
species climatic envelope, that climatic optima for growth
differ among local populations, correlating with their
climate of origin, and that climatic tolerance only partially
overlap between *provenances*. These tests also suggest
that populations located at the extremities of the natural
distribution inhabit climates that are suboptimal for their
growth and development. For population 2 (coming from
latitude 59.1°N), as illustrated in A, the optimal climate
(2.5°C) is warmer than the climate of its geographic source
(-2.5°C), while for southern populations the opposite
pattern occurs. The overall picture is a negative correlation
between Δ and the latitude of origin of the population as
shown in panel B by the example in the case of *Pinus
contorta* (Rehfeldt *et al.* 1999). This pattern is consistent
with theoretical predictions that asymmetric gene flow
from the core to peripheral populations increases
maladaptation at the edges of the natural distribution
(Kirkpatrick & Barton, 1997; Garcia-Ramos & Kirkpatrick
1997, see also text).

<figure><img src="image_68.png" /><figcaption>Annual temperature at testing site (°C)</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Latitude of origin of the populations</figcaption></figure>
TABLE 1

<table><caption>Examples of observed LD pollen and seed dispersal in trees (more than 3 km for pollen and 1 km for seeds). The Table is arranged first by propagule type (pollen or seed), than by vector type (wind, insects, birds, bats, elephants, fish or generally vertebrates) and dispersal type (potential, viable or effective, defined in a footnote below), and finally alphabetically by species name.</caption><thead><tr><th rowspan="2">Species</th><th colspan="3">Dispersal system</th><th rowspan="2">Location</th><th rowspan="2">Method</th><th colspan="2">Dispersal Distance</th><th rowspan="2">Reference</th></tr><tr><th>Propagule</th><th>Vector</th><th>Type<sup>a</sup></th><th>Maximum</th><th>Proportion ≥ threshold<sup>b</sup></th></tr></thead><tbody><tr><td>Betula spp.</td><td>Pollen</td><td>Wind</td><td>Potential</td><td>Central-North-Eastern Europe</td><td>Aerobiologic and phenological analysis</td><td>1000 km</td><td></td><td>Siljamo et al. 2008</td></tr><tr><td>Pinus banksiana and Picea glauca</td><td>Pollen</td><td>Wind</td><td>Potential</td><td>Canada</td><td>Aerobiologic analysis</td><td>3000 km</td><td></td><td>Campbell et al. 1999</td></tr><tr><td>Pinus sylvestris</td><td>Pollen</td><td>Wind</td><td>Viable</td><td>Northern Europe</td><td>Aerobiologic and phenological analysis</td><td>600 km</td><td></td><td>Varis et al. 2009</td></tr><tr><td>Pinus taeda</td><td>Pollen</td><td>Wind</td><td>Viable</td><td>Eastern North America</td><td>Aerobiologic analysis</td><td>40 km</td><td></td><td>Williams 2010</td></tr><tr><td>Cecropia obtusifolia</td><td>Pollen</td><td>Wind</td><td>Effective</td><td>Central America</td><td>Genetic paternity analysis</td><td></td><td>10% > 14 km</td><td>Kaufman et al. 1998</td></tr><tr><td>Fraxinus excelsior</td><td>Pollen</td><td>Wind</td><td>Effective</td><td>Scotland</td><td>Genetic parentage analysis</td><td></td><td>25-35% > 3 km</td><td>Bacles et al. 2006</td></tr><tr><td>Pinus sylvestris</td><td>Pollen</td><td>Wind</td><td>Effective</td><td>Spain</td><td>Genetic mixture analysis</td><td></td><td>4.3% > -100 km</td><td>Robledo-Arnuncio 2011</td></tr><tr><td>Quercus robur</td><td>Pollen</td><td>Wind</td><td>Effective</td><td>Eastern Europe</td><td>Genetic parentage analysis</td><td></td><td>35% > 80 km</td><td>Buschbom et al. 2011</td></tr><tr><td>Populus trichocarpa</td><td>Pollen</td><td>Wind</td><td>Effective</td><td>Western North America</td><td>Genetic paternity analysis</td><td></td><td>5% > -5-10 km</td><td>Slavov et al. 2009</td></tr><tr><td>Ficus spp.</td><td>Pollen</td><td>Insects</td><td>Effective</td><td>Central America</td><td>Genetic parental reconstruction</td><td>14 km (isolated mother trees)</td><td></td><td>Nason et al. 1998</td></tr><tr><td>Ficus sycomorus</td><td>Pollen</td><td>Insects</td><td>Effective</td><td>Namibia</td><td>Genetic paternity analysis</td><td>165 km</td><td></td><td>Ahmed et al. 2009</td></tr><tr><td>Sorbus domestica</td><td>Pollen</td><td>Insects</td><td>Effective</td><td>Central Europe</td><td>Genetic paternity analysis</td><td></td><td>-1% > 12-16 km</td><td>Kamm et al. 2009</td></tr><tr><td>Swietenia humilis</td><td>Pollen</td><td>Insects</td><td>Effective</td><td>Central America</td><td>Genetic paternity analysis</td><td></td><td>40-80% ≥ 4 km (in small fragments)</td><td>White et al. 2002</td></tr><tr><td>Fraxinus excelsior</td><td>Seed</td><td>Wind</td><td>Effective</td><td>Scotland</td><td>Genetic parentage analysis</td><td>1.4 km</td><td>46-53% > 3 km</td><td>Bacles et al. 2006</td></tr><tr><td>Annona glabra</td><td>Seed</td><td>Birds</td><td>Potential</td><td>Australia</td><td>Empirically-based simulations of vector movements and seed passage time</td><td>5.2 km</td><td>1% > 4 km</td><td>Westcott et al. 2008</td></tr><tr><td>Xylopia hypolampra and 7 other species</td><td>Seed</td><td>Birds</td><td>Potential</td><td>Cameroon</td><td>Empirically-based simulations of vector movements and seed passage time</td><td>6.9 km</td><td></td><td>Holbrook & Smith 2000</td></tr><tr><td>Ficus carica and Morus alba</td><td>Seed</td><td>Bats</td><td>Potential</td><td>Israel</td><td>Empirically-based simulations of vector movements and seed passage time</td><td>20 km</td><td>17% > 1 km</td><td>Tsoar et al. 2011</td></tr><tr><td>Tamarindus indica</td><td>Seed</td><td>Elephants</td><td>Potential</td><td>Myanmar (Burma)</td><td>Empirically-based simulations of vector movements and seed passage time</td><td>5.4 km</td><td>50% > 1.2 km</td><td>Campos-Arceiz et al. 2008</td></tr><tr><td>Duroia duckei and 2 other species</td><td>Seed</td><td>Fish</td><td>Potential</td><td>Peru</td><td>Empirically-based simulations of vector movements and seed passage time</td><td>5.5 km</td><td>5% > 1.7 km</td><td>Anderson et al. 2011</td></tr><tr><td>Prunus mahaleb</td><td>Seed</td><td>Vertebrates</td><td>Potential</td><td>Spain</td><td>Genetic maternal analysis</td><td></td><td>33% >1500 m</td><td>Jordano et al. 2007</td></tr><tr><td>Sorbus domestica</td><td>Seed</td><td>Vertebrates</td><td>Effective</td><td>Central Europe</td><td>Genetic paternity analysis</td><td>12.2 km</td><td></td><td>Kamm et al. 2009</td></tr></tbody></table>

<footnote><sup>a</sup>Three types are distinguished: "potential" dispersal is the distance dispersed by a propagule (pollen or seeds) at any, commonly unknown, condition; "viable" is the same as "potential" but excluding non-viable propagules; cases of "effective" dispersal are the pollen that gave rise to seeds, or seeds that established, yielding seedlings, saplings or young/adult plants.</footnote>

<footnote><sup>b</sup>The proportion (in %) of propagules dispersed to equal or greater distances than the specified threshold. The threshold distances were defined by the authors of each study, often arbitrarily or according to features of the study landscape and/or populations.</footnote>
effects of gene flow on adaptation in trees. We conclude
that the positive effects of gene flow may often dominate
its negative effects, although regional variation may
influence the balance of those effects. We finally elaborate
on the theoretical and experimental approaches that
should be implemented to improve our ability to predict
the scale and distribution of gene flow effects on forest
ecosystems in the context of climate change.

How far do seeds and pollen disperse
in trees?

Gene flow in plants is mediated by both seed and pollen
dispersal, which vary greatly among species (Ennos 1994).
Seed and pollen dispersal however have distinct effects on
the rate of demographic spread, and the rate at which
genes move across the range of a species. The spatial scale
of effective propagule dispersal (glossary) in trees
depends on a variety of physical and biological processes

that determine the amount and availability of pollen and
seeds, their movement, their viability before and during
movement, and the probability of successful pollination
leading to viable seed and the seedling establishment
rates. Different combinations of these components may
yield *effective dispersal* (glossary) distances spanning
from a few centimeters to thousands of kilometers
(Nathan *et al.* 2008), generally following markedly
leptokurtic patterns. Aerobiological studies show that
airborne tree pollen (both viable and non-viable) has the
potential to be transported in substantial amounts over
hundreds to thousands of kilometers (Table 1, Fig. 2).
However, documented dispersal distances of viable pollen
(yet prior to successful fertilization) are about one order
of magnitude shorter, up to 600 km (Table 1). Documented
distances of effective pollen dispersal (when pollination
led to successful mating) are of lower magnitude, up to
100 km (Table 1). Documented wind-driven effective seed
dispersal is up to a few kilometers (Table 1), thus about two

<caption>FIGURE 2</caption>

Virtual long distance pollen dispersal of *Pinus taeda*

<figure><img src="image_70.png" /></figure>

Virtual pollen release, using the Regional Atmospheric Modeling System (RAMS) and its Eulerian-Lagrangian particle transport module (HYPACT). This regional atmospheric simulation was forced with meteorological data from National Oceanic and Atmospheric Administration (NOAA) NCEP-DOE Reanalysis II data set. The experimental settings are described in Bohrerova et al. (2009). The figure shows a portion of the southeast United States, centered on eastern South Carolina. Pollen was arbitrarily released from two locations, in North Carolina outer banks (black point) and South Carolina (grey point), at a simulated afternoon on 27 March 2006, corresponding

with the peak of pollen release at the Duke forest, NC.
The dispersing pollen plumes (black for NC pollen,
grey for SC pollen) are shown as "clouds." The wind
was moderate, mainly toward the northeast. The figure
shows a snapshot of the pollen plume at 6:00 AM, 36
hours after the release. The viability is resolved by the
model as an additional property of the pollen.
Mortality due to UV and vapour pressure deficit is
calculated, with rates fitted to empirical equations,
based on observations of a bench-scale experiment.
The pollen in the image, 36 hours after release is about
40% viable.
orders of magnitude shorter than effective pollen dispersal. Although animal-mediated seed dispersal can reach a scale of tens of kilometers, pollen dispersal distances are in general considerably longer than that of seeds, and especially in wind-driven dispersal systems.

Detecting mating events over hundreds of kilometers is experimentally very difficult and many biological processes take place (viability, phenology, pollen competition) between pollen release and fertilization. Similarly, quantifying the multiple demographic and dispersal parameters affecting seed dispersal and subsequent establishment is still hard to accomplish (but see Nathan et al. 2011). These experimental constraints and biological reality both explain the disparities between reported potential and effective dispersal distances. Their relative contributions remain unresolved and addressing this point will require innovation (see section 5.1). Interactions between physical and biological processes influencing long-distance dispersal (LDD) and the way these will be affected in a changing climate are particularly poorly understood. The local environment around the release point of the propagule and the conditions at which the propagule was released can have strong effects of either decreasing or increasing the rate and extent of LD propagule transport (Bohrer et al. 2008; Wright et al.

2008). Also important is the interaction between weather patterns that drive extreme LDD, and the conditions that affect survival during the dispersal event, particularly in pollen, where a migration event may span many hours and days (Schueler et al. 2005).

In that respect, predicting how climate change itself may modify patterns of gene flow in the future remains a formidable challenge. Projected changes such as warmer atmosphere (Kuparinen et al. 2009), or higher fecundity, earlier maturation, and changes (reduction or increase) in surface wind speed (Nathan et al. 2011) can be implemented into mechanistic models of plant spread in future environments. Although this approach facilitates assessing the upper and lower bounds of future gene flow and population spread, it is unlikely to provide accurate predictions for a particular species and system, due to the large uncertainty about key hard-to-measure parameters, such as the spatial patterns of phenological schedules for pollen and of post-dispersal survival of seeds. Moreover, genetic variation in traits affecting dispersal of both pollen and seeds is common in plants and, accordingly, plant dispersal traits have been observed to evolve fast in response to environmental change, especially in the context of range expansion (Darling et al. 2008). In particular, increased frequency of traits facilitating seed

<caption>FIGURE 3</caption>

<figure><img src="image_71.png" /><figcaption>Effect of dispersal distance on the evolutionary load</figcaption></figure>

The average distance between population mean and optimal phenotype (i.e. maladaptation) is predicted from equation 2.2 in Box 2, as $\sqrt{2 \cdot V_s \cdot L}$. The model is illustrated with the example of the evolution of *bud set* date in the Sitka spruce as explained in Box 2. The standard deviation within populations for phenotypic expression of bud-set date is 17 days. *Heritability* was taken to be 0.5. No data is available on the strength of

*stabilizing selection* $V_s$. Compiling the results of many studies (e.g. as shown in Johnson & Barton 2005), the median value of $P/V_s$ is 0.2 (solid line). We also used a stronger value of *stabilizing selection* ($P/V_s$=1, dashed line). The optimal dispersal distance is smaller when selection is stronger (and thus the spatial fitness gradient steeper) but the *evolutionary load* is also much reduced (compare dashed and solid line).
dispersal have been found in recently founded populations, at the expanding edge of the range, while the converse trend was found in fragmented southern populations at the rear end (Riba et al. 2009). Identifying how global change (including climate change) will alter selection pressures acting on dispersal is therefore crucial in order to predict the extent of gene flow in future environments.

Despite gaps that prevent us from precisely predicting the extent of LD gene flow, the available data on effective dispersal and predicted habitat shifts still suggest that these two processes may operate over comparable scales for many tree species, but also that inter-specific variation

in the rate and magnitude of LD gene flow may substantially affect the variation in the response of forest trees to climate change (Nathan et al. 2011, Table 1).

## What are the possible effects of gene flow on adaptation?

We can conceptualize the potential effects of gene flow on adaptation through a simple model of species range evolution, where selection varies both in space and time (Pease et al. 1989; Polechova et al. 2009; Box 2). Local climate can be thought of as imposing specific selection pressures on a complex set of phenotypic traits (e.g.

## BOX 2

### A simple conceptual framework illustrating the antagonistic effects of gene flow

Two closely related theoretical models (Pease et al. 1989; Polechova et al. 2009) have explored the question of evolution in environments changing both in space and time, mimicking the effects of climate change in species with wide distributions. These models envision a species distributed along some linear environmental gradient, such as the Sitka spruce distributed along a large latitudinal gradient of temperatures (Mimura & Aitken 2007a). We use this empirical example to illustrate potentially realistic values of parameters for the model. We assume fitness depends quadratically on how well an individual tree is adapted to its local ecological conditions, i.e., how its phenotype matches the local optimum. By averaging over phenotypes in the local population, one may then write the mean fitness (here the Malthusian population exponential growth rate) in a given location as:

$$ \bar{r} = r_0 - \frac{P}{2V_s} - \frac{(\bar{z} - \theta)^2}{2V_s}, \quad (2.1) $$

where $r_0$ is the contribution to population growth of an individual with the optimal phenotype, $V_s$ describes how well individuals that deviate from this optimal phenotype perform (and is thus inversely related to the strength of stabilizing selection), $\theta$ is the optimal phenotype in that location, $\bar{z}$ is the local mean phenotype in the population and $P$ is the local phenotypic variance around this mean. This expression shows that the mean fitness in a variable population subject to stabilizing selection is reduced in two ways: (1) Standing load ($P/2V_s$): caused by phenotypic variation and present even when the mean phenotype matches the optimum; (2) Evolutionary load [$\bar{z} - \theta]^2/2V_s$]: caused by departure of the mean phenotype from the local optimum (Lande & Shannon 1996). Evolutionary loads can be generated by selection that varies over space or time (see Bridle et al. 2009 for a review).

Pease et al. (1989) and Polechova et al. (2009) make specific predictions about how migration might affect the evolutionary load in a changing climate. Their models assume the optimum phenotype changes linearly through space, with slope $b$. This is similar to Sitka spruce, where bud set date increases linearly with the local mean annual temperature, which itself varies linearly with the distance to the Southern margin of the species range: assuming current bud set date corresponds to the optimum, this

gives an estimate of $b$ such that optimal bud set date increases by 13 days every °C, or by 3.24 days per 100 km (Mimura & Aitken 2007a; Aitken et al. 2008). Climate change can be approximated by this gradient of optimal phenotypes being constantly shifted through space at rate $v$. According to different climate models, mean annual temperature may increase by 3-5°C in the generation time of Sitka Spruce (Aitken et al. 2008), which gives an estimate of $v$ as a shift of approximately 1000 to 2000 km per generation. Migration is modelled as a diffusion process, with average distance between parent and offspring $\sigma$. There are feedbacks between the evolution of the mean phenotype through time and space and that of the population density, mediated through gene flow and the local growth rates (Pease et al. 1989; Polechova et al. 2009).

Further assuming that genetic variation for the trait under selection is relatively weak and does not vary through space, Pease et al. (1989) predict that the loss of fitness at the scale of the range due to evolutionary load is approximately

$$ L \approx \frac{1}{2} \left( \frac{\sigma b}{\sqrt{V_s}} + \frac{v^2}{\sigma^2} - \frac{G}{V_s} \right) \quad (2.2) $$

where $G$ is the genetic variance for the trait. Though this prediction might be crude in the case of forest trees with large within population genetic variance, it has heuristic value. Indeed, the first term in parentheses can be interpreted as the component of phenotypic mismatch due to spatial variability in the optimal phenotype and gene flow (migration load); this part of the load increases with dispersal distance ($\sigma$). The second term describes phenotypic mismatch due to the lagging response of the mean phenotype to temporal change in the local optimum (lag load). This part of the load decreases with dispersal distance because migration helps the species track its shifting optimum through space. The third term in (2.2) shows that the evolutionary load declines with genetic variance $G$ because response to selection increases. Dispersal distance also affects the evolution of genetic variance $G$ (Barton 2001; Polechova et al. 2009; Bridle et al. 2010), with positive effects on the evolutionary load (equation 2.2), but negative effects on the standing load $P/2V_s$ (equation 2.1)
phenology, frost hardiness, growth, seed size), and
defining different optimal trait values through various
trade-offs, depending on specific combinations of climatic
conditions encountered within the range. As an illustration,
in Sitka Spruce, trees originating from higher latitude with
lower annual mean temperature cease growing earlier in
the season than trees from lower latitude, when grown in
common garden (Mimura & Aitken 2007a), suggesting
different optimal *bud set* (glossary) dates along
temperature gradients within the range. The simple
conceptual model in Box 2 connects adaptation to
demography by assuming that an individual's contribution
to population growth declines as it departs from the locally
optimal phenotype. This model suggests that gene flow
has antagonistic effects on adaptation by modifying the
various sources of *genetic load* (glossary) depressing
population mean fitness, and thus population growth (see
Bridle *et al.* 2009 for a review). We here review these
effects by considering their alternative evolutionary
consequences.

Gene flow constrains local adaptation

Because gene flow homogenises allele frequencies across
space, high gene flow could constrain *adaptive divergence*
(glossary) along environmental gradients (Garcia-Ramos
& Kirkpatrick 1997; Bohrer *et al.* 2005; but see Barton 2001;
Yeaman & Guillaume 2009; Bridle *et al.* 2010 for a revised
consideration of the strength of such constraints). Some
theoretical models predict in particular that gene flow
from large central populations into small peripheral ones
may swamp local adaptation in marginal areas, preventing
range spread beyond some critical environmental limit
(Kirkpatrick & Barton 1997; review in Bridle & Vines 2007).
Gene flow then causes *phenotypic clines* (glossary) for
adaptive traits to deviate from their optima (Box 2).
Interestingly, the interaction of strong gene flow with
selection on multiple traits could result in some *phenotypic
clines* being flatter, and some steeper, than optimal
because of genetic and selective interactions among traits
(Guillaume 2011, Duputié *et al.* in revision). Furthermore,
the constraining effects of migration on divergence are
predicted to be more severe when divergence involves
many loci of small effects rather than few major genes with
large effects on the phenotype (Yeaman & Guillaume
2009). In forest trees, the former situation seems to be the
most common (Neale & Kremer 2011).

Common garden experiments in forest trees (Box 1)
suggest that genotypes can perform poorly when
transferred to climates far from their location of origin.
Maladaptation of LD migrants could thus reduce the mean
fitness in forest tree populations, generating a *migration
load* (glossary). Such *migration load* would be of concern
if gene flow is extensive over long distances (see section
2) and if phenotypic mismatch of immigrants is mostly due
to long lasting genetic effects (Aitken et al. 2008). In a
Swedish population of *Pinus sylvestris*, Nilsson (1995)

indeed found that offspring sired by naturally dispersing
pollen had significantly slower growth and higher freezing
resistance than expected if offspring were sired only by
pollen produced locally. Natural pollination thus resulted
in a phenotypic shift corresponding to that expected if
most pollen originated from higher latitudes by 1 to 2
degrees (Nilsson 1995).

There is however little evidence that gene flow has
strongly limited adaptation in forest trees in the past.
Comparison of genetic differentiation at neutral molecular
markers versus *adaptive traits* (glossary) repeatedly
suggests that extensive gene flow (presumably mostly
through pollen dispersal) has not prevented rapid adaptive
divergence of extant populations (Savolainen *et al.* 2007;
Kremer *et al.* 2010, for theoretical predictions see Kremer
& Le Corre 2011). Populations under different climates may
however have diverged while still being far from the locally
optimum phenotype. Determining how much observed
*phenotypic clines* deviate from what would be optimal
under local conditions is however difficult to assess, and
remains an open issue in evolutionary biology (Barton,
2001; Butlin *et al.* 2003). When *provenance tests* (glossary)
have been replicated over a broad range of climatic
conditions, *provenance* (glossary) *reaction norms*
(glossary, Box 1) allow comparison of optimal and original
climate for each population. Mismatches are not
uncommon (Wang *et al.* 2010) and interestingly some
studies (Rehfeldt *et al.* 1999; Rehfeldt *et al.* 2002) found
more of them at the edge of distributions: e.g. populations
of *Pinus contorta* from locations with extreme climates
grow better in milder conditions, closer to the core, than
in their original location. Such a pattern is consistent with
the theoretical expectation that gene flow from the core
increases maladaptation in marginal populations (Garcia-
Ramos & Kirkpatrick 1997).

Gene flow enhances the response to selection

Natural selection operates by sieving from genetic
variation found within populations. Local genetic diversity
is therefore the fuel of evolutionary change. Forest tree
populations harbour high diversity both for molecular
markers and quantitative traits (Hamrick *et al.* 1992), with
*heritabilities* (glossary) typically above 0.4 for wood
characteristics or phenological traits such as bud set date
(Cornelius, 1994). The maintenance of such high levels of
quantitative variation for traits closely linked to fitness
remains a paradox where strong *stabilizing selection*
(glossary) is acting to reduce variation within populations
(see Johnson & Barton 2005 for a review). Theoretical
models predict that increases in *genetic variance*
(glossary) due to gene flow could be substantial (Barton,
2001). Along climatic gradients, the increase in genetic
variance within localities due to gene flow is predicted to
be proportional to the change in mean *breeding value*
(glossary) along the typical dispersal distance. For
example, in Sitka Spruce, the *breeding value* for *bud set*
date varies by 3.24 days every 100 km while the within-
population phenotypic standard deviation of *bud set* date
is typically about 10-25 days (Mimura & Aitken 2007a,
Aitken et al. 2008); LD gene flow over distances of about
100 km would then lead to *heritability* for *bud set* date
greater than 0.4 even for relatively strong *stabilizing*
*selection*. In addition to the mean dispersal distance, the
whole shape of the *dispersal kernel* (glossary) is predicted
to affect the spatial distribution of genetic variation and
diversity of migrants (Travis et al. 2010; Fayard et al. 2009),
especially in the context of range expansion.

Extensive gene flow in trees is generally thought of as a
major explanation for their high within-population diversity
(Hamrick et al. 1992). Together with the strong selection
acting at the juvenile stage experienced in trees, this may
allow rapid adaptation to changing climate without large
significant reductions in population mean fitness (for
empirical examples of rapid genetic changes in forest trees
see Jump et al. 2006). There is however little direct
empirical demonstration of this. Using a mechanistic
model of beech stand dynamics, Kramer et al. (2008)
predicted little effect of pollen dispersal distance on the
evolution of within-stand genetic diversity, but their model
ignored the potentially large phenotypic divergence of
immigrants (e.g. Nilsson 1995). If gene flow between

differentiated populations is a persistent source of *genetic
variation*, there should be strong correlations between
genetic variance within populations and the amount of
heterogeneity in the environment at the regional scale. In
a study of 142 populations of lodgepole pine, Yeaman &
Jarvis (2006) indeed found that the variance for growth
among individuals from the same locality measured in
*provenance tests* (glossary, Box 1) correlated with regional
climatic heterogeneity.

Gene flow facilitates the tracking
of environmental change

Assuming an initially locally adapted population, the new selection pressures induced by climate change will cause the optimal phenotype to deviate from the extant mean phenotype, thus creating a *lag load* (glossary, Box 2). Gene flow will reduce the *lag load* in two different ways: (i) by helping track the shift of the optimum phenotype through dispersal of pre-adapted genotypes found somewhere else in the range (see equation 2.2 in Box 2), and (ii) by augmenting the *response to selection* (glossary) through a general increase in genetic variation (see previous paragraph). Quite generally, dispersal should help adaptation in environments that are changing both in space and time (Blanquart & Gandon 2011). Beyond

<caption>FIGURE 4</caption>

Long distance effective pollen dispersal of *Pinus sylvestris* L.

<figure><img src="image_72.png" /></figure>

Estimated effective pollen immigration rates into a
*Pinus sylvestris* remnant (encircled with dashed line)
from five long-distant populations (encircled with
continuous lines) in Central Spain, obtained using
maximum-likelihood genetic mixture analysis
combined with Monte Carlo assessment of small

parameter uncertainty. Continuous (resp. dashed)
arrows indicate pollen immigration rates significantly
LTE
(resp. non-significantly) different from zero. 95%
confidence intervals between brackets (modified from
Robledo-Arnuncio 2011 with permission from The
Publisher).
theoretical predictions, there is little experimental evidence
on the impact of gene flow on the *lag load* of tree
populations engaged in an evolutionary race with a
changing environment. Simulating the evolution of growth
cessation date in pine and birch, Kuparinen *et al.* (2010)
found that the large lag that developed after one century
of climate warming (about 40 days behind the optimal
date) was only moderately reduced (by 2 days) by the
higher pollen and seed dispersal distances of birch along
the climatic gradient.

Gene flow affects inbreeding levels

Gene flow may also affect genetic variation for fitness
independently from the issue of adaptation to climate, by
affecting genetic resemblance between mates. Early
*inbreeding depression* (glossary) is widespread in largely
outcrossing species such as trees and inbreeding was
found more frequently in isolated marginal populations
(Mimura & Aitken 2007b), which may depress their mean
fitness and their ability to persist in a changing
environment. In many plant species with small isolated
populations, progeny formed by crosses between
populations show higher fitness than that obtained by
within population crosses (e.g. Willi & Fischer 2005). Such
a pattern of *heterosis* (glossary) is however not expected
to be generally very strong in trees, due to the
combination of large population size and extensive gene
flow (for theoretical predictions see Glémin *et al.* 2003;
Lopez *et al.* 2009). Conversely, small amounts of gene flow
between formerly isolated populations can also in theory
seriously disrupt fitness due to negative interactions
between genes having evolved separately (Edmands &
Timmerman 2003), but contrasting results from artificial
crosses between distant populations fail to provide solid
evidence for this type of *outbreeding depression*
(glossary) in trees (Harfouche *et al.* 2000; Goto *et al.* 2011).
Because tree populations are seldom isolated from each
other, *outbreeding depression* due to negative gene
interactions is expected to be rare (Frankham *et al.* 2011).

How will the different effects of gene
flow balance each other in the context
of climate change?

We here suggest that the positive effects of gene flow
may often dominate negative effects for forest trees
confronted to climate change. This is in particular the case
due to the specific life cycle of forest trees. However, this
balance of effects is likely to be modulated by (i) the
regional context (e.g. expanding edge or retracting part of
the range), and (ii) the dispersal syndrome (relative
strength of pollen versus seed dispersal).

Balance between antagonistic effects

Maladaptation in a changing climate is caused by
mismatches between optimal and realized mean
phenotypes, due to environments that are changing in

time and space too fast for the population to adjust to
these changes. Such mismatch depresses the mean fitness
of populations, generating an *evolutionary load* (glossary,
Box 2). Setting aside the effects of gene flow on *genetic
variance*, the simple model in Box 2 predicts that there is
an optimal level of migration that minimizes such
*evolutionary load* under a shifting climate: when dispersal
distance is short, the *lag load* decreases fast with
increasing migration, which helps the population track the
optimal climate (Fig. 3). If gene flow is too high, however,
local adaptation is prevented (*migration load*) and
maladaptation increases (albeit slowly) with increasing
migration (Fig. 3). The optimal dispersal distance is higher
if the environment changes more quickly in time in a given
location, and if selection varies less sharply in space (Box
2). As an illustration, in the evolution of bud set date in
Sitka Spruce (Mimura & Aitken 2007a; Aitken et al. 2008,
see Box 2), the optimal migration distance is relatively
large (immigrants should on average originate from
locations with mean temperature differing by more than
3°C to the local site, i.e. more than 1000km). This suggests
that, for a range of realistic dispersal distances, the positive
tracking effect of dispersal should dominate its negative
effects on local adaptation.

Once the effects of gene flow on the evolution of *genetic variance* are taken into account, the constraining effects of migration on adaptation in marginal populations is much weakened: very high gene flow seems instead to facilitate adaptation across a wide array of environmental conditions, but at the cost of a reduced fitness everywhere in the range, which could ultimately compromise species persistence (for theoretical predictions see Barton 2001; Polechova *et al.* 2009), and this effect on population fitness becomes greater when the stochastic effects of finite populations are included (Bridle *et al.* 2010). The very high fecundity, long life span and strong competition at the juvenile stage, which are characteristic of forest trees, could in principle permit very high *genetic load* causing massive mortality at the juvenile stage, without having much impact on adult density. Further exploration of connections between forest trees population dynamics and genetic diversity are needed to conclude when the demographic cost of adaptation compromises persistence.

Overall, models that integrate different antagonistic
effects predict that intermediate levels of gene flow (e.g.
between one and ten migrants per generation) suffice to
replenish *genetic variance* eroded by drift and selection,
and alleviate *inbreeding depression* without causing large
*migration load*, thus maximize mean fitness in
heterogeneous environments (Lopez *et al.* 2009; see also
Blanquart & Gandon 2011 for the case of spatio-temporal
variation). Empirical evidence that migration enhances
fitness in marginal habitats of several plant species
supports such predictions (Kawecki 2008). We lack similar
direct evidence in forest trees. The admixture of genotypes
of diverse geographical origin is increasingly thought of as