Sara Branco
email: sara.mayer.branco@gmail.com

Ke Bi

Hui-Ling Liao

Pierre Gladieux

Helene Badouin

Christopher E Ellison

Nhu H Nguyen

Rytas Vilgalys

Kabir G Peay

John W Taylor

Hélène Badouin

Thomas D Bruns

Hélène Badouin

Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes

Keywords: Adaptation, population genomics, mycorrhizal fungi, Suillus brevipes

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The environment in which organisms live plays an important role in shaping the distribution of genetic variability within and between populations across time and space. Environmental adaptation in natural populations leaves genomic signatures of positive selection that shed light on the origins and maintenance of genetic variation within and between populations [START_REF] Dobzhansky | Genetics of Natural Populations. XVI. Altitudinal and Seasonal Changes Produced by Natural Selection in Certain Populations of Drosophila Pseudoobscura and Drosophila Persimilis[END_REF][START_REF] Novembre | Spatial patterns of variation due to natural selection in humans[END_REF]. Natural selection affects the levels of genetic variability and can lead to selective sweeps, where sites linked to a selected mutation show reduced within-population variability [START_REF] Nielsen | Molecular Signatures of Natural Selection[END_REF]. One way to studying environmental adaptation is to detect genomic regions under positive selection. This approach can be done by using large-scale single nucleotide polymorphism (SNP) datasets, detecting statistical outliers and characterizing genes under selection. Another method to studying environmental adaptation is ecological association, in which one detects allele frequencies that exhibit significant statistical association with environmental variables that are used as proxies for natural selection [START_REF] Frichot | Testing for Associations between Loci and Environmental Gradients Using Latent Factor Mixed Models[END_REF][START_REF] Frichot | LEA: An R package for landscape and ecological association studies[END_REF]. Genome scans for environmental adaptation produce lists of genetic polymorphisms potentially involved in local adaptation that can be compared to genomic regions under positive selection and provide a better assessment of environmental adaptation.

Studies on the genomics of environmental adaptation are very often conducted in biological systems with pre-defined ecotypes that are recognized based on obvious morphological traits, ecological preferences, or both (e.g. [START_REF] Jones | The genomic basis of adaptive evolution in threespine sticklebacks[END_REF], Laurent et al. 2015[START_REF] Nachman | The genetic basis of adaptive melanism in pocket mice[END_REF][START_REF] Schweizer | Genetic subdivision and candidate genes under selection in North American grey wolves[END_REF]. Such studies aim to document the genomic basis of ecotypic differences and are based on a priori ideas on relevant population units and environmental parameters that might cause genetic Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

divergence. For most taxa, including the vast majority of microbiological diversity, a priori information on ecotypic variation is scarce if available at all, making studies of environmental adaptation more challenging to implement. Fungi in particular offer little obvious phenotypic variation potentially linked to fitness differences at the population level and there are only a few documented cases of fungal local adaptation relating to environmental differences, such as heavy metal contaminated soils (Colpaert et al. 2011), genetically different plant pathogenic hosts (e.g. [START_REF] Salvaudon | Genetic diversity in natural populations: a fundamental component of plant-microbe interactions[END_REF][START_REF] Enjalbert | Genetic evidence of local adaptation of wheat yellow rust (Puccinia striiformis f. sp. tritici) within France[END_REF]) and thermal clines [START_REF] Mboup | Genetic structure and local adaptation of European wheat yellow rust populations: the role of temperature-specific adaptation[END_REF]. Recent ecological genomic approaches have revealed additional evidence of fungal adaptation related to subtle environmental variation and allowed candidate genomic regions involved in environmental adaptation to be identified. Such regions provide indirect evidence for fungal adaptation to environmental conditions (an approach coined 'reverse ecology'; [START_REF] Li | Reverse Ecology" and the Power of Population Genomics[END_REF] and can be the basis for developing tests aimed at revealing the mechanisms responsible for the process of fungal adaptation. Relevant examples include ecotypic variation in the saprobe Neurospora crassa linked to temperature gradients [START_REF] Ellison | Population genomics and local adaptation in wild isolates of a model microbial eukaryote[END_REF], adaptation to cheese environment in the genus Penicillium [START_REF] Ropars | Adaptive Horizontal Gene Transfers between Multiple Cheese-Associated Fungi[END_REF], host specificity in the wheat pathogen Zymoseptoria graminicola (Stukenbrock et al. 2011) and adaptation to soil chemistry in the ectomycorrhizal fungus Suillus brevipes (Branco et al. 2015). The latter study focused on the same species that we examine here and documented two recently diverged and isolated populations in an obligate fungal symbiont of pine from coastal and montane California. The two populations showed very limited genetic differentiation, but a gene homologous to Nha-1, a membrane Na + /H + exchanger involved in salt tolerance in plants and

Saccharomyces cerevisiae [START_REF] Nass | Intracellular Sequestration of Sodium by a Novel Na + /H + Exchanger in Yeast Is Enhanced by Mutations in the Plasma Membrane H + -ATPase. Insights into mechanisms of sodium tolerance[END_REF][START_REF] Apse | Salt Tolerance Conferred by Overexpression of a Vacuolar Na+/H+ Antiport in Arabidopsis[END_REF]) differed between the two populations and showed evidence of positive selection. This report was the first to find salt as a potential driver for ectomycorrhizal fungal population differentiation, even though soil chemistry is long known to affect fungal physiology and growth [START_REF] Smith | Mycorrhizal Symbiosis[END_REF].

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017

Accepted Article

This article is protected by copyright. All rights reserved.

Here we investigated the genomic basis of adaptation to climate at the continental scale in S. brevipes by expanding the genomic dataset studied of Branco et al. (2015). Adding to the previously documented coastal and montane California populations, we sequenced the genomes of individuals from five additional sites in Canada (two locations were sampled), Colorado, Minnesota and Wyoming (United States of America), totaling 55 S. brevipes individuals. This expanded sample encompassed a wide range of environments and enabled us to test whether S. brevipes genomes display signatures of adaptation to distinct climatic regimes. Specifically, we asked whether there is evidence for local adaptation in this species across North America by detecting selective sweeps across the S. brevipes genome and assessing whether large-scale climatic parameters such as temperature and precipitation are associated with genetic differentiation. We hypothesized that this wind-dispersed species would show both strong population structure across the continent and a pattern of climatic adaptation, with genomic signatures of positive selection in regions correlated to temperature and/or precipitation and involved in environmental response. We found strong population differentiation across North America as well as significant genotype-climate associations in genomic regions under strong positive selection. These regions included genes involved in transmembrane transport and helicase activity that are known from other systems to be involve with cold stress response.

Materials and methods

Isolate collection, culturing and DNA extraction and sequencing

We collected fruitbodies of Suillus brevipes and cultured dikaryotic individuals from them from seven pine forests selected from the western and northern range of the species. Specifically we sampled two sites in California (USA), two sites in Alberta (Canada), and one site, each, in Wyoming, Colorado and Minnesota (USA, Fig. 1). A total of 55 dikaryotic S. brevipes individuals were included in our study. This sample contained 28 dikaryotic S. brevipes cultures from California that we had been previously sequenced and reported on Branco et al. (2015), and 27 new individuals from the other

Accepted Article

This article is protected by copyright. All rights reserved.

parts of range. We sequenced the whole genomes of the new cultures with pair-end 100bp Illumina technology following the protocol described in Branco et al. (2015). We combined these new data with the sequences generated for the S. brevipes populations from coastal and montane California studied in Branco et al. (2015) to produce a total sample of 55 S. brevipes individuals. We used four individuals as outgroups, one S. luteus individual and three individuals of Suillus sp. (S106, S108, S110), a species closely related to S. pseudobrevipes (Nguyen et al, in press). Suillus luteus reads were obtained from the reference genome (http://genome.jgi.doe.gov/Suilu1/Suilu1.home.html) and we sequenced the Suillus sp. (S106, S108, S110) whole genomes from fruitbodies collected in Florida and processed in the same manner as S. brevipes individuals. All details on collection sites for all individuals are summarized in Table S1. Raw reads were deposited in Short Read Archive and accession numbers can also be found in Table S1.

Read filtering, mapping and variant calling

We used a Galaxy workflow to perform read filtering and mapping [START_REF] Goecks | Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences[END_REF], Blankenberg et al. 2010[START_REF] Giardine | Galaxy: a platform for interactive large-scale genome analysis[END_REF]; Supplement 1). In short, we filtered the reads using Trim Galore!

Accepted Article

This article is protected by copyright. All rights reserved.

frequencies and population genetic parameters. ANGSD takes base calling, mapping and alignment errors into account, making it ideal to analyze the genomes included in this study. This software, which was developed for low and medium coverage data and diploid organisms, is equally useful for studying dikaryons.

Population structure and individuals genealogy

We used RAxML (Stamatakis 2014) to infer the S. brevipes individuals maximum-likelihood genealogy based on the full set of SNPs (using the model GTRGAMMA and 100 bootstrap replicates). Population structure was studied using principal components analysis (PCA) with ngsCovar implemented in ngsTools [START_REF] Fumagalli | Quantifying Population Genetic Differentiation from Next-Generation Sequencing Data[END_REF] and NgsAdmix [START_REF] Skotte | Estimating individual admixture proportions from next generation sequencing data[END_REF], a software for assigning individuals to differentiated genetic clusters and inferring admixture of related individuals using genotype likelihood data. Analyses were performed with the number of clusters, K, ranging from two to seven and with five replicates per K value. Average likelihoods were compiled using the Evanno method [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] to identify genetically homogeneous groups of individuals. The divergence metrics F ST and D xy were computed for the homogeneous populations only (the two California populations, Colorado and Minnesota) using ANGSD and ngsStat, respectively.

Summary statistics and linkage disequilibrium

Summary statistics and linkage disequilibrium decay were estimated for the four homogenous populations (the two California populations, Colorado and Minnesota). We computed nucleotide diversity (π), Watterson theta (θ) and Tajima's D with ANGSD [START_REF] Korneliussen | ANGSD: Analysis of Next Generation Sequencing Data[END_REF] and heterozygosity for all populations using ngsStat in ngsTools [START_REF] Fumagalli | Quantifying Population Genetic Differentiation from Next-Generation Sequencing Data[END_REF]. To assess linkage disequilibrium (LD) decay, we calculated the correlation coefficient (r2) between any two loci across

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

the genome using genotypes. r2 was calculated using VCFTools 0.1.13 [START_REF] Danecek | The variant call format and VCFtools[END_REF]) with a maximum of 5000bp between SNPs.

Association between genotypes and environmental variables

Genotype-environment associations were assessed using Latent Factor Mixed Models 1.4 (LFMM, [START_REF] Frichot | Testing for Associations between Loci and Environmental Gradients Using Latent Factor Mixed Models[END_REF]), a univariate mixed-model approach that estimates the effects of hidden factors representing background residual levels of population structure. We acquired environmental parameters for each S. brevipes homogeneous population using geo-referenced environmental layer datasets [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]

Accepted Article

This article is protected by copyright. All rights reserved.

In LFMM the optimal latent factor number (K) influences the false discovery rate and should be selected based on the mean genomic inflation factor (λ), which is the median of the squared z scores divided by 0.456. Based on 5 LFMM runs, K=6 had the λ estimates closer to 1.0 for both envPC1 and envPC2. We obtained the median z scores of all five runs and re-adjusted the p-values using the Benjamini-Hochberg correction procedure and set the false discovery rate to 10% as recommended in the program manual. We selected the top 5% SNPs significantly correlated with envPC1 and envPC2

to investigate gene ontology enrichment (see below).

Genome scans for selective sweeps

We used SweeD 3.2.11 [START_REF] Pavlidis | SweeD: Likelihood-based detection of selective sweeps in thousands of genomes[END_REF] to search for selective sweeps. These genomic regions are defined by reduced genetic variation (with an excess of rare variants and of high frequency derived variants) around mutation sites due to strong and recent positive selection that has led to the fixation of a single allele [START_REF] Nielsen | Molecular Signatures of Natural Selection[END_REF]. SweeD implements a composite likelihood ratio (CLR) test based on the SweepFinder algorithm [START_REF] Nielsen | Genomic scans for selective sweeps using SNP data[END_REF]. The CLR uses the variation of the whole or derived site frequency spectrum of a genomic sequence to compute the ratio of the likelihood of a selective sweep at a given position to the likelihood of a null model without selective sweep. The null hypothesis relies on the site frequency spectrum (SFS) of the whole genomic sequence rather than on a standard neutral model, which makes it more robust to demographic events such as population expansions [START_REF] Pavlidis | SweeD: Likelihood-based detection of selective sweeps in thousands of genomes[END_REF], Nielsen et al. 2005). SweeD was run both separately for each of the populations identified using analyses of population structure to detect population-specific selective sweeps, and on all samples to detect species-wide selective sweeps. We used the empirical derived site frequency spectra of the whole genome as a background for all scaffolds, using the four outgroup individuals to orient SNPs. We considered only scaffolds where the first and last SNPs were more than 3 kb distant, and computed CLR every 1 kb.

Accepted Article

This article is protected by copyright. All rights reserved.

CLR values not belonging to the 95th percentile were considered as outliers. As the CLR is not designed to test if close outlier values belong to the same selective sweep, we simply defined genomic intervals affected by selective sweeps by merging consecutive outlier positions into single intervals and adding 500 bp at the flanks of each interval.

Outlier Gene Ontology (GO) enrichment analysis

We conducted gene ontology (GO) enrichment analyses on selective sweeps and LFMM outliers using ClueGO 2.2.4, a Cytoscape plugin Java tool that extracts knowledge on putative function for large clusters of genes. We modified the default settings to include minimal number and percentage of genes for GO term/pathway selection and the minimal and maximal level for GO Tree interval as one and twelve respectively. Significance of GO terms was calculated with a two-sided hypergeometric test and p-values corrected using the Benjamini-Hochberg approach. We also identified enrichment of functional domains associated to the two environmental principal components. We used a hypergeometric test to assess differences in functional domains associated with envPC1 and envPC2.

The geometric test was further applied to each domain with genetic variation regarding envPC1 and envPC2, with a total of 862 domains examined. Individual functional domains with FDR higher or equal than 0.05 were likely observed by chance after correcting for multiple comparisons.

Results

Suillus brevipes genomes

We obtained over one billion raw Illumina reads and over 900 million high-quality reads from 55 S.

brevipes individual genomes, that yielded 1,495,484 high quality single nucleotide polymorphisms across the 52.7Mb genome. Polymorphism level in this species was close to 3% segregating sites.

Table S1 summarizes all information on per-individual number of generated raw reads, number of high-quality reads, alignment rate and genome coverage (which varied between 21x and 110x).

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison

Accepted Article

This article is protected by copyright. All rights reserved.

Population structure, individual genealogy and summary statistics

We found distinct and clearly differentiated S. brevipes populations across North America, including the two Californian populations previously described in Branco et al. (2015). The genealogy of individuals based on the full set of SNPs (Fig. 2) showed the Canadian samples as paraphyletic and sister to all other populations, with Minnesota nested within Castle Rock and Whitecourt basal to

Wyoming, Colorado and coastal and montane California.

The ngsAdmix analyses were consistent across distinct K values, revealing coastal California, montane California, Colorado and Minnesota as genetically homogeneous (Fig. 3, Figure S2).

Individuals from the two populations in Canada and Wyoming showed substantial membership in clusters identified in other populations suggesting potential admixture (Fig. 3).

The SNP-based PCA also showed S. brevipes population genetic structure (Fig. 4), with the two Californian populations previously analyzed appearing as clearly differentiated but accounting for a small portion of the genomic variability encompassed in the sampled populations. Apart from Minnesota and Colorado, the remaining populations were not as clearly differentiated. Individuals from the two Canadian sites and Wyoming showed some degree of overlap in the ordination space and Canadian individuals spread over both PC1 and PC2 indicating high genetic variability.

California, Colorado and Minnesota populations showed similar per-base nucleotide diversity (π),

Watterson theta (θ) and heterozygosity (Table 1). Linkage disequilibrium (LD) decay curves were very similar across populations (Fig. 5), showing a strong sharp decay in the first hundreds of base pairs. Such pattern indicates very small linkage blocks and therefore pervasive outcrossing in this species.

Genomic pairwise F ST revealed high levels of differentiation between populations of S. brevipes (Table 2). The California montane and coastal populations were the least differentiated (with F ST CAc- CAm = 0.09), while Minnesota was the most differentiated, with F ST values of 0.3 and over when compared to California and Colorado. These results are in accordance with both the phylogeny and admixture analyses.

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

Selective sweep detection and association between genomic sites and environmental variables

Selective sweep analyses revealed over 1500 genomic regions with signatures of positive selection in each population and all populations combined. Only 95 sites were recovered as swept in all analyses.

Sweeps were prevalent across the majority of scaffolds, with larger scaffolds containing more swept sites (Fig. S3). Gene ontology (GO) enrichment analyses failed to detect enriched GO categories across the sweeps, indicating that sites under positive selection span a wide range of functions.

LFMM analyses found 18744 sites significantly associated with envPC1 (which is heavily influenced by winter temperature and precipitation) and 6776 sites with envPC2 (which is correlated with weather heterogeneity throughout the year). GO enrichment analyses on these sites revealed that the categories significantly enriched for envPC1 were ATPase-driven transmembrane transporters and ATPase-dependent helicases (Fig. 6a), while envPC2 was significantly enriched with protein transmembrane transporter activity (Fig 6b).

We found a total 862 functional domains in the LFMM sites correlated with the two environmental principal components (727 in envPC1, 459 in envPC2, with an overlap of 324 domains). A hypergeometric test revealed significant differences between the domains in envPC1 and envPC2 (p = 1.8e-93). Of 862 functional domains with genetic variation either in envPC1 or envPC2, three were identified as not randomly enriched in envPC1 or/and envPC2 (Table 3; FDR < 0.05; p < 0.05). The domains' functions were protein kinase, WD40 repeats and nucleic acid binding and DNA/RNA helicase.

We found 43 individual genomic sites that were both significantly associated with environmental parameters and showed signatures of positive selection as assessed by LFMM and selective sweep analysis (Table 4). These sites were located across the scaffolds and most were outliers restricted to a single population. These results suggest the existence of local adaptation, with different alleles being selected for in different populations. Eleven genomic sites currently have no predicted function and four sites were included in start/stop codons. The latter were adjacent to transporter genes and hydrolases, indicating potential loss-of-function mutations. Notably, the Na + /H + exchanger associated

Accepted Article

This article is protected by copyright. All rights reserved.

with salt tolerance found to be fixed in the coastal California population in Branco et al. (2015) was also found to be a relevant genomic region for that population only, and it was under strong positive selection and strongly associated with both envPC1 and envPC2.

Discussion

Suillus brevipes biogeography

Suillus brevipes is an obligate mycorrhizal associate of Pinus subgenus Pinus and is widespread in North America (Smith, A. H. & Thiers, H. D. 1964, Nguyen et al., in press), tightly associated with the range of its hosts. This species is one of the most ruderal of the genus and is commonly found with young pine trees at forest margins and also with introduced and invasive pines (Ashkannejhad & Horton 2006[START_REF] Peay | A strong species-area relationship for eukaryotic soil microbes: island size matters for ectomycorrhizal fungi[END_REF][START_REF] Hynson | Identities and distributions of the coinvading ectomycorrhizal fungal symbionts of exotic pines in the Hawaiian Islands[END_REF]). This high level of host specificity limits the species distribution and creates a mosaic continental structure of habitat for S. brevipes.

In spite of the prolific production of aerially dispersed meiotic spores in Suillus spp. [START_REF] Peay | Measuring ectomycorrhizal fungal dispersal: macroecological patterns driven by microscopic propagules[END_REF][START_REF] Dahlberg | Size, distribution and biomass of genets in populations of Suillus bovinus (L.: Fr.) Roussel revealed by somatic incompatibility[END_REF], there is mounting evidence that effective dispersal is limited to a fairly fine spatial scale, which is consistent with the idea that 'almost nothing is everywhere' (Gladieux et al. 2015a, Taylor et al. 2006[START_REF] Peay | Not every fungus is everywhere: scaling to the biogeography of fungal-plant interactions across roots, shoots and ecosystems[END_REF]. In fact, recent work showed that the spore loads of ectomycorrhizal fungi, including Suillus species, decrease exponentially as spores travel away from their source, support the hypothesis that these fungi are not very effective at dispersing over long distances [START_REF] Peay | Measuring ectomycorrhizal fungal dispersal: macroecological patterns driven by microscopic propagules[END_REF]. For example, even though S. pungens was estimated to produce 8 trillion spores/km 2 within its host forest, its spore load decreased exponentially with distance from the forest edge, with spore levels insufficient to consistently colonize pine seedlings at distances greater than 1km [START_REF] Peay | Measuring ectomycorrhizal fungal dispersal: macroecological patterns driven by microscopic propagules[END_REF]. Such observed dispersal limitation at the local scale is also consistent with both the population genomics results of Branco et al. (2015) revealing two isolated and locally adapted S. brevipes populations in California and the results reported here. Long distance spore dispersal in S. brevipes is thus likely to be rare, and if it occurs, immigrant genotypes may be less fit

Accepted Article

This article is protected by copyright. All rights reserved.

than local ones. This explanation is consistent with our finding that North American S. brevipes is subdivided into multiple, clearly circumscribed genetic clusters.

The least distinct and most genetically diverse populations occur in the two Canadian sites and to a less extent in Wyoming. The finding of individuals assigned to multiple clusters in Canada might suggest admixture between individuals from differentiated populations. This hybridization idea is consistent with the known postglacial history of P. banksiana and P. contorta, the only hosts in this area. In fact, these pine species form their own hybrid zone in western Alberta within the Whitecourt site that we sampled (Yeatman 1967). Population genetic studies [START_REF] Godbout | Phylogeographic structure of jack pine (Pinus banksiana; Pinaceae) supports the existence of a coastal glacial refugium in northeastern North America[END_REF]) and pollen records (Williams et al. 2004) suggest that these pine species met in Alberta as the glaciers retreated, by reinvasion from Eastern North America (P. banksiana) and from one or more western mountain corridors (P. contorta) (Critchfield 1985).

Alternatively, Canadian admixture can be interpreted as Alberta hosting the oldest S. brevipes populations, and therefore the most diverse, from which all other populations would be derived. This interpretation is consistent with the individuals genealogy and with previous findings of artifactual inference of admixture in ancestral populations in other organisms [START_REF] Kauer | Nonneutral admixture of immigrant genotypes in African Drosophila melanogaster populations from Zimbabwe[END_REF][START_REF] Schöfl | Microsatellite variation and differentiation in African and non-African populations of Drosophila simulans[END_REF][START_REF] Gladieux | On the Origin and Spread of the Scab Disease of Apple: Out of Central Asia[END_REF]), but would require the presence of a pine refugia in this area during glaciation..

Adaptation in Suillus brevipes

We found evidence consistent with abiotic environmental adaptation in S. brevipes populations across North America. Genes found both under positive selection and significantly associated with climatic environmental variables were mainly involved in transmembrane transport of substances and helicase activity. Even though these gene categories are involved in multiple cell functions, we argue that they can be directly related to environmental adaptation, namely in the physiological responses to conditions determined by variations in temperature and precipitation.

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

Transmembrane transporter genes are known to play important roles, including responding to environmental stress in wide variety of organisms (Chan et al. 2011[START_REF] Krulwich | Molecular aspects of bacterial pH sensing and homeostasis[END_REF][START_REF] Mansour | The plasma membrane transport systems and adaptation to salinity[END_REF], Chen et al. 2003[START_REF] Apse | Salt Tolerance Conferred by Overexpression of a Vacuolar Na+/H+ Antiport in Arabidopsis[END_REF][START_REF] Rizzello | Protein cold adaptation strategy via a unique sevenamino acid domain in the icefish (Chionodraco hamatus) PEPT1 transporter[END_REF][START_REF] Mahajan | Cold, salinity and drought stresses: an overview[END_REF][START_REF] Heide T Van Der | Osmoregulated ABC-transport system of Lactococcus lactis senses water stress via changes in the physical state of the membrane[END_REF]. In fungi, transporters are involved in environmental sensing and are crucial for detecting and responding to a wide variety of cues, including chemistry, temperature and water availability [START_REF] Aguilera | Cold response in Saccharomyces cerevisiae: new functions for old mechanisms[END_REF], Bahn & Mühlschlegel 2006, Bahn et al. 2007[START_REF] Maldonado-Mendoza | A Phosphate Transporter Gene from the Extra-Radical Mycelium of an Arbuscular Mycorrhizal Fungus Glomus intraradices Is Regulated in Response to Phosphate in the Environment[END_REF][START_REF] Javelle | High-affinity ammonium transporters and nitrogen sensing in mycorrhizas[END_REF][START_REF] Harrison | A Phosphate Transporter from Medicago truncatula Involved in the Acquisition of Phosphate Released by Arbuscular Mycorrhizal Fungi[END_REF][START_REF] Dietz | The aquaporin gene family of the ectomycorrhizal fungus Laccaria bicolor: lessons for symbiotic functions[END_REF][START_REF] Pettersson | Aquaporins in yeasts and filamentous fungi[END_REF][START_REF] Peter | Ectomycorrhizal ecology is imprinted in the genome of the dominant symbiotic fungus Cenococcum geophilum[END_REF]. Temperature and precipitation differences across the sampled sites can therefore be directly linked to the enriched transporter gene categories detected in our study. Cold stress in particular affects the activity of membrane transporters due to changes in membrane fluidity that can have multiple nefarious consequences, including decreased nutrient uptake [START_REF] Aguilera | Cold response in Saccharomyces cerevisiae: new functions for old mechanisms[END_REF]. Water stress on the other hand can be mitigated by membrane proteins that mediate water transfer (aquaporins, Carbrey & Agre 2009). These transporters are known to occur in fungi [START_REF] Pettersson | Aquaporins in yeasts and filamentous fungi[END_REF], being involved in the transport of both water and ammonia in ectomycorrhizal fungi and presumed to be directly linked to mechanisms of environmental adaptation [START_REF] Dietz | The aquaporin gene family of the ectomycorrhizal fungus Laccaria bicolor: lessons for symbiotic functions[END_REF][START_REF] Peter | Ectomycorrhizal ecology is imprinted in the genome of the dominant symbiotic fungus Cenococcum geophilum[END_REF].

Helicases are motor proteins that separate nucleic acids strands. These enzymes are involved in many fundamental processes such as replication, transcription, translation and nucleic acid repair and are known to be involved in environmental stress response. RNA helicases play a role in responding to abiotic stressors by changing the export of nuclear mRNA, translation initiation, mRNA decay, rRNA processing, transcription and the progression of cell cycle [START_REF] Owttrim | RNA helicases and abiotic stress[END_REF]. RNA helicases are known to be involved in adaptation to cold in bacteria and fungi (Barria et al. 2013[START_REF] Pandiani | Differential involvement of the five RNA helicases in adaptation of Bacillus cereus ATCC 14579 to low growth temperatures[END_REF][START_REF] Aguilera | Cold response in Saccharomyces cerevisiae: new functions for old mechanisms[END_REF]). In the aforementioned "reverse ecology" study of Neurospora crassa populations [START_REF] Ellison | Population genomics and local adaptation in wild isolates of a model microbial eukaryote[END_REF]), an RNA helicase was one of several genes hypothesized to be important to cold adaptation and one of only two genes whose role in cold adaptation could be confirmed by gene deletion. Cold stress can lead to the stabilization of RNA secondary structure, impairing transcription, translation and degradation, and RNA helicases act to melt such structures

Accepted Article

This article is protected by copyright. All rights reserved.

allowing for RNA degradation (Barria et al. 2013). In Saccharomyces cerevisiae, RNA helicases are induced during the early stages of cold stress response [START_REF] Aguilera | Cold response in Saccharomyces cerevisiae: new functions for old mechanisms[END_REF].

DNA helicases are involved in DNA repair, as they allow repair machinery to access the single stranded DNA [START_REF] Goldman | The DNA Damage Response in Filamentous Fungi[END_REF], and are also induced by cold stress (Vashisht et al. 2005), allowing for repair mechanisms to act on this molecule. Further functional tests addressing the role of candidate genes are needed to investigate potential differences in environmental adaptation (including cold and water stress response) in S. brevipes individuals across North America.

The membrane Na + /H + exchanger described in Branco et al. (2015) to be under strong positive selection, was again found here among the highly differentiated genes and strongly associated with envPC1 and envPC2 in the coastal California population. This gene is homologous to Nha-1, known for enhancing salt tolerance in plants and S. cerevisiae [START_REF] Apse | Salt Tolerance Conferred by Overexpression of a Vacuolar Na+/H+ Antiport in Arabidopsis[END_REF][START_REF] Nass | Intracellular Sequestration of Sodium by a Novel Na + /H + Exchanger in Yeast Is Enhanced by Mutations in the Plasma Membrane H + -ATPase. Insights into mechanisms of sodium tolerance[END_REF]) and a good candidate for local adaptation in S. brevipes. It is therefore not surprising that this gene is located in a putative selective sweep region and located in a D xy outlier genomic window. However, our LFMM analyses were restricted to climatic variables and did not include soil parameters, so possible explanations for the strong correlation with envPCs include pleiotropic functional properties or linkage to another gene truly correlated with weather variability.

Concluding remarks

Our study provides a unique large-scale view of within species evolution and sheds light on the factors driving environmental adaptation in fungi. We used a large-scale dataset to detect genomic regions correlated to climatic parameters and harboring signatures of positive selection without any prior information on phenotypic differences, population units or environmental factors that could underlie potential adaptation. Such reverse ecology approach is very suitable for microbial systems, where there is very little a priori information. Fungi in particular have small genome sizes, simple morphologies, short generation times and a multitude of life styles, being exceptional models for studying adaptation in eukaryotes [START_REF] Gladieux | Fungal evolutionary genomics provides insight into the mechanisms of adaptive divergence in eukaryotes[END_REF]. Reverse ecology has in fact revealed hidden

Accepted Article

This article is protected by copyright. All rights reserved.

fungal environmental adaptation in several species so far (e.g. Branco et al. 2015[START_REF] Ellison | Population genomics and local adaptation in wild isolates of a model microbial eukaryote[END_REF], Gladieux et. al. 2015b and Badouin et. al. submitted to this issue). Genome scans are therefore a useful framework for identifying putative adaptive genes responsible for local adaptation, and the results can motivate functional tests aimed at linking genotypes to phenotypes and are crucially needed for a full understanding of the genetic basis of adaptation.

Figures

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. Only the genes that were significantly correlated with one envPC compared to the other were listed. 'Percentage of genes in GO term' refers to the percentage of genes from the uploaded cluster that were associated with the term, compared to all S. brevipes genes associated with the same term. The numbers on bars correspond to the numbers of genes from the analyzed cluster associated with each category. Gray -helicase activity; black -transmembrane transporters; white -other categories.

  , C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Figure 1 -

 1 Figure 1 -Sampling localities of the studied Suillus brevipes populations in North America. Redcoastal California (16 ind); purple -montane California (9 individuals); yellow -Canada (Alberta), Whitecourt (5 individuals); orange -Canada (Alberta), Castle Rock (7 individuals); blue -Wyoming (7 individuals); black -Colorado (5 individuals); green -Minnesota (6 individuals).

Figure 3 -

 3 Figure 3 -Estimated S. brevipes population structure as inferred by NGSAdmix. Each individual is represented by a bar partitioned into K color segments (K=2 to K=7). CAc -coastal California; CAmmontane California; CANc -Canada (Alberta), Castle Rock; CANw -Canada (Alberta), Whitecourt; CO -Colorado; MN -Minnesota; WY -Wyoming.

Figure 5 -

 5 Figure 5 -Linkage disequilibrium decay in each population across the whole genome. Populations are color-coded as in Fig. 1: Red -coastal California; purple -montane California; black -Colorado; green -Minnesota.

  ., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Figure 6 -

 6 Figure6-Percentage and number of genes assigned to different gene ontology categories significantly correlated with envPC1 (a) and envPC2 (b) as assessed by ClueGO (p<0.005). Only the genes that were significantly correlated with one envPC compared to the other were listed. 'Percentage of genes in GO term' refers to the percentage of genes from the uploaded cluster that were associated with the term, compared to all S. brevipes genes associated with the same term. The numbers on bars correspond to the numbers of genes from the analyzed cluster associated with each category. Gray -helicase activity; black -transmembrane transporters; white -other categories.

  transmembrane transporter activity Hydrolase activity, acting on acid anhydrides, catalyzing transmembrane movement of substances Ion transmembrane transporter activity Primary active transmembrane transporter activity Macromolecule transmembrane transporter activity Protein transmembrane transporter activity P-P-bond-hydrolysis-driven transmembrane transporter activity P-P-bond-hydrolysis-driven protein transmembrane transporter activity ATPase activity, coupled to transmembrane movement of substances Protein-transmembrane transporting ATPase activity ATPase coupled ion transmembrane transporter activity ATPase activity, coupled to movement of substances Single-stranded DNA-dependent ATP-dependent DNA helicase

  Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

  Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

  Figure 2 -Rooted maximum-likelihood genealogy inferred from the full SNP dataset. We used one individual of S. luteus (Slu) and three individuals of Suillus sp. as outgroups (in gray).

	Colorado
		Wyoming
		Montane
		California
		Coastal
		California
		Castle Rock
	Whitecourt	(Canada)
	(Canada)
	Minnesota
		The remaining
	colors correspond to Fig. 1: Red -coastal California; purple -montane California; yellow -Canada
	(Alberta), Whitecourt; orange -Canada (Alberta), Castle Rock; blue -Wyoming; black -Colorado;
	green -Minnesota. * indicates bootstrap support value >90.

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (

2017

). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Acknowledgements

We thank Cathy Cripps and her lab for providing Suillus cultures and Matthew Smith for inadvertently providing an extra outgroup species. Holly Edes and Ravi Alla provided technical assistance, Lucinda Lawson assisted in obtaining and analyzing climatic data, Olivier François helped implementing LFFM, Steven Wu assisted in implementing enrichment analyses. Tatiana Giraud provided comments on earlier drafts and enthusiastically supported this study. Three anonymous reviewers provided insightful comments on this manuscript. This work was funded by NSF Dimension of Biodiversity grant DBI-1046115 and DEB-1257528, and used the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley, supported by NIH S10 Instrumentation Grants S10RR029668 and S10RR027303.

Accepted Article

This article is protected by copyright. All rights reserved. Ashkannejhad S, Horton TR (2006) Ectomycorrhizal ecology under primary succession on coastal sand dunes: interactions involving Pinus contorta, suilloid fungi and deer. The New Phytologist, 169, 345-354. Badouin H, Gladieux P, Gouzy J, et al. Submitted. Widespread selective sweeps throughout the genome of the model pathogenic fungi and identification of effector candidates. Molecular Ecology. Bahn Y-S, Mühlschlegel FA (2006) CO2 sensing in fungi and beyond. Current Opinion in Microbiology, 9, 572-578. Bahn Y-S, Xue C, Idnurm A et al. (2007) Sensing the environment: lessons from fungi. Nature Reviews Microbiology, 5, 57-69. Barria C, Malecki M, Arraiano CM (2013) Bacterial adaptation to cold. Microbiology, 159, 2437-2443. Blankenberg D, Von Kuster G, Coraor N et al. (2010) Galaxy: a web-based genome analysis tool for experimentalists. Current Protocols in Molecular Biology, FM Ausubel (ed), pp 1-21. Wiley Branco S, Gladieux P, Ellison CE et al. (2015) Genetic isolation between two recently diverged populations of a symbiotic fungus. Molecular Ecology, 24, 2747-2758. Carbrey JM, Agre P (2009) Discovery of the Aquaporins and Development of the Field. In: Aquaporins Handbook of Experimental Pharmacology, PDE Beitz (ed), pp. 3-28. Springer Berlin Heidelberg. Chan CX, Reyes-Prieto A, Bhattacharya D (2011) Red and Green Algal Origin of Diatom Membrane Transporters: Insights into Environmental Adaptation and Cell Evolution. PLOS ONE, 6, e29138. Chen D, Toone WM, Mata J et al. (2003) Global Transcriptional Responses of Fission Yeast to Environmental Stress. Molecular Biology of the Cell, 14, 214-229. Version postprint Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. (2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved. Colpaert JV, Wevers JHL, Krznaric E, Adriaensen K (2011) How metal-tolerant ecotypes of ectomycorrhizal fungi protect plants from heavy metal pollution. Annals of Forest Science, 68, 17-24. Critchfield WB (1985) The late Quaternary history of lodgepole and jack pines. Canadian Journal of Forest Research, 15, 749-772.

Version postprint

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Version postprint

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

Data accessibility

DNA sequences: Short Read Archive accessions in Table S1 (Supporting information).

Version postprint

Comment citer ce document : Branco, S., Bi, K., Liao, H.-L., Gladieux, P., Badouin, H., Ellison, C. E., Nguyen, N. H., Vilgalys, R., Peay, K. G., Taylor, J. W., Bruns, T. D. ( 2017). Continental-level population differentiation and environmental adaptation in the mushroom Suillus brevipes. Molecular Ecology, 26 (7), 2063-2076. DOI : 10.1111/mec.13892

Accepted Article

This article is protected by copyright. All rights reserved.

Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.p5dm3

Supporting information

Table S1 -Suillus brevipes isolate location, pine host, coordinates and elevation, number of raw reads, average genome coverage, and Short Read Archive accession numbers.

Table S2 -Loadings for each principal component (envPC) in the environmental PCA for the studied Suillus brevipes populations.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.