# Population genomics of a forest fungus reveals high gene flow and climate adaptation signatures

Jørn Henrik Sønstebø, Emiliano Trucchi, Jenni Nordén, Inger Skrede, Otto Miettinen, Sajeet Haridas, Jasmyn Pangilinan, Igor Grigoriev, Francis Martin, Håvard Kauserud, et al.

► To cite this version:

Jørn Henrik Sønstebø, Emiliano Trucchi, Jenni Nordén, Inger Skrede, Otto Miettinen, et al. Population genomics of a forest fungus reveals high gene flow and climate adaptation signatures. Molecular Ecology, 2022, 31 (7), pp.1963-1979. 10.1111/mec.16369. hal-04001302

HAL Id: hal-04001302

https://hal.inrae.fr/hal-04001302

Submitted on 22 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
ORIGINAL ARTICLE

MOLECULAR ECOLOGY

WILEY

# Population genomics of a forest fungus reveals high gene flow and climate adaptation signatures

Jørn Henrik Sønstebø<fnref n="1" /> ID

Emiliano Trucchi<fnref n="2" /> ID

Jenni Nordén<fnref n="3" /> ID

Inger Skrede<fnref n="1" /> ID

Otto Miettinen<fnref n="4" /> ID

Sajeet Haridas<fnref n="5" /> ID

Jasmyn Pangilinan<fnref n="5" /> ID

Igor V. Grigoriev<fnref n="5" />,<fnref n="6" /> ID

Francis Martin<fnref n="7" /> ID

Håvard Kauserud<fnref n="1" /> ID

Sundy Maurice<fnref n="1" /> ID

**Correspondence**

Jørn Henrik Sønstebø and Sundy Maurice, Section for Genetics and Evolutionary Biology (EVOGENE), Department of Biosciences, University of Oslo, Oslo, Norway.

Emails: Jorn.H.Sonstebo@usn.no (J.H.S.); sundymaurice@gmail.com (S.M.)

**Funding information**

Norges Forskningsråd, Grant/Award Number: 254746

**Abstract**

Genome sequencing of spatially distributed individuals sheds light on how evolution structures genetic variation. Populations of *Phellopilus nigrolimitatus*, a red-listed wood-inhabiting fungus associated with old-growth coniferous forests, have decreased in size over the last century due to a loss of suitable habitats. We assessed the population genetic structure and investigated local adaptation in *P. nigrolimitatus*, by establishing a reference genome and genotyping 327 individuals sampled from 24 locations in Northern Europe by RAD sequencing. We revealed a shallow population genetic structure, indicating large historical population sizes and high levels of gene flow. Despite this weak substructuring, two genetic groups were recognized; a western group distributed mostly in Norway and an eastern group covering most of Finland, Poland and Russia. This substructuring may reflect coimmigration with the main host, Norway spruce (*Picea abies*), into Northern Europe after the last ice age. We found evidence of low levels of genetic diversity in southwestern Finland, which has a long history of intensive forestry and urbanization. Numerous loci were significantly associated with one or more environmental factors, indicating adaptation to specific environments. These loci clustered into two groups with different associations with temperature and precipitation. Overall, our findings indicate that the current population genetic structure of *P. nigrolimitatus* results from a combination of gene flow, genetic drift and selection. The acquisition of similar knowledge especially over broad geographic scales, linking signatures of adaptive genetic variation to evolutionary processes and environmental variation, for other fungal species will undoubtedly be useful for assessment of the combined effects of habitat fragmentation and climate change on fungi strongly bound to old-growth forests.

**KEYWORDS**

demography, forest fragmentation, genotype-environment associations, local adaptation, *Picea abies*, post-glacial colonization, wood decay fungi

<footnote n="1">¹Section for Genetics and Evolutionary Biology (EVOGENE), Department of Biosciences, University of Oslo, Oslo, Norway</footnote>

<footnote n="2">²Department of Life and Environmental Sciences, Marche Polytechnic University, Ancona, Italy</footnote>

<footnote n="3">³Norwegian Institute for Nature Research, Oslo, Norway</footnote>

<footnote n="4">⁴Finnish Museum of Natural History, University of Helsinki, Finland</footnote>

<footnote n="5">⁵US Department of Energy Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA</footnote>

<footnote n="6">⁶Department of Plant and Microbial Biology, University of California Berkeley, Berkeley, California, USA</footnote>

<footnote n="7">⁷Université de Lorraine, INRAE, UMR Interactions Arbres/Microorganismes, INRAE GrandEst-Nancy, Champenoux, France</footnote>
# 1 INTRODUCTION

The evolutionary mechanisms shaping the genetic variation of any species are determined by the biology of the species and the environment in which it lives. In Northern Europe, genetic variation was highly influenced by the last period of glaciation, particularly for species that survived this period in refugia, subsequently colonizing their current habitats via different colonization routes (Hewitt, 1999). During the last period of glaciation, many species retreated to refugia in Southern or Eastern Europe, often with small effective population sizes, resulting in low levels of genetic variation within refugia, but large differences between refugia. After recolonization of the areas previously covered by ice, genetic variation was in many species highly influenced by gene flow from different source refugia (Hewitt, 1996). Genetic variation within the distribution range of a species can also be impacted by the spatial variation in selection pressure, leading to local adaptation and high levels of genetic differentiation between populations (Savolainen et al., 2013). The level of differentiation depends on the strength and type of selection, and the amount of gene flow between populations (Tigano & Friesen, 2016).

Human activities have become a major driving force shaping the genetic variation of natural populations, leading to fragmentation of large populations, locally decreases in genetic diversity and greater differentiation between populations (Young et al., 1996). Decreases in genetic diversity may reduce both the potential for migration to new areas and the ability to adapt to a changing climate. Thus, understanding how species respond to changes in environmental conditions is of the utmost importance, in both theoretical and applied ecological genomics, to facilitate adaptive management practices.

Wood-inhabiting fungi are pivotal for the functioning of the forest ecosystems, due to their ability to decompose organic material. The occurrence of wood-inhabiting fungal species is driven by local factors relating to the amount and quality of dead wood substrate available, as well as by regional factors relating to landscape history, old-growth forest connectivity and climate (Abrego et al., 2017; Nordén et al., 2013, 2018; Penttilä et al., 2006). It remains largely unknown how these different factors interact to shape the intraspecific genetic variation and population structure of wood-inhabiting fungi. Previous population genetic studies of wood-inhabiting fungi have suggested that population sizes and regional gene flow levels are both high, resulting in weak genetic differentiation (e.g Högberg et al., 1999; Kauserud & Schumacher, 2002). Adaptation potential is often higher in populations with large effective population sizes, but, high levels of gene flow between populations inhabiting areas with different climates would be expected to reduce the ability to adapt (Tigano & Friesen, 2016). Decreasing area and connectivity of forests rich in dead wood due to land-use change and modern forestry practices have resulted in smaller, spatially distant local populations for many wood-inhabiting fungi. Fragmentation may limit gene flow and favor genetic drift due to small population sizes (Young & Clarke, 2000), decreasing genetic diversity and increasing genetic differentiation between populations, as reported for some

rare wood-inhabiting fungi (Högberg & Stenlid, 1999; Penttilä et al., 2006).

We studied the wood-inhabiting fungus *Phellopilus nigrolimitatus* in Northern Europe to determine how different evolutionary processes have shaped the genetic diversity of this red-listed species. *Phellopilus nigrolimitatus* is an outcrossing basidiomycete (Hymenochaetaceae) found predominantly in old-growth forests, where it fruits on large logs in relatively late (3 and 4) decay stages (Kauserud & Schumacher, 2001). *Phellopilus nigrolimitatus* is one of the many fungi affected by modern forestry practices in Fennoscandia, including clear-cuttings, resulting in a decrease in population sizes in recent decades (Henriksen & Hilmo, 2015; Stokland & Kauserud, 2004). By contrast, in areas with a more continuous habitat, the genetic structure of populations probably results from the postglacial expansion of the main substrate of this fungus, Norway spruce, in Northern Europe. From eastern glacial refugia in the Russian plains, Norway spruce expanded westward in Northern Europe along two main routes: a northern route extending from Russia into Finland and central Sweden and Norway, and a southern route extending from Russia into the Baltic states, southern Sweden and Norway, stopping along the mountains separating western and eastern Norway (Chen et al., 2019; Tollefsrud et al., 2008, 2009). The main range expansion of Norway spruce reached its current western limit relatively recently (about 2000–3000 years ago, Tollefsrud et al., 2008).

The distribution of *P. nigrolimitatus* spans a large range of temperature and precipitation, from a continental climate in Russia to an oceanic climate in Western Europe. These regional differences in climate may result in different selection pressures, with the potential to generate differences in allele frequencies at genetic loci affected by selection. Conversely, as *P. nigrolimitatus* seems to have reached its current distribution relatively recently, there may not have been sufficient time for selection to produce strong patterns of local adaptation. However, in many studies focusing on introduced species, differences in adaptation in spatially separated populations have been shown to evolve within only a few generations (e.g van Boheemen et al., 2019; Koskinen et al., 2002; Medley et al., 2019).

We explored the population genetic processes that have shaped the genome-wide variation in *P. nigrolimitatus*, by analysing restriction-site associated DNA (RAD) sequencing data for 24 Northern Europe populations, with main emphasis on populations in Norway and Finland. We hypothesized that the genetic structure of *P. nigrolimitatus* results principally from the range expansion of its main host, *P. abies*, after the last glaciation period. We also expected to find a decline in genetic diversity and an increase in genetic differentiation between locations in areas of extensive forestry and human land use. Given the relatively large environmental differences in our study area, we hypothesized that selection would have generated significant genotype-environment associations in *P. nigrolimitatus*. As shown in previous studies of climate adaptation in fungi (Branco et al., 2017; Dal Grande et al., 2017; Ellison et al., 2011), we expected temperature and precipitation to be important selective pressures driving the distribution of genetic variation at relevant genomic loci in *P. nigrolimitatus*.
2 | MATERIALS AND METHODS

2.1 | Study organism and sampling

The basidiomycete *Phellopilus nigrolimitatus* (family Hymenochaetaceae) is known to have a heterothallic sexual system, implying that mating occurs between haploid mycelia of different mating types (Kauserud & Schumacher, 2001). Sexual reproduction occurs in the fruiting structure, the sporocarp. After a short diploid (2n) phase, meiosis generates haploid basidiospores (n) that germinate into homokaryotic mycelia. When two compatible homokaryotic mycelia from the same or different sporocarps mate, a fertile dikaryotic mycelium (n + n) is formed, in which nuclear exchange and reciprocal migration of the nuclei occur (Boddy, 2016). This dikaryotic mycelium corresponds to the dominant stage of an individual and, in this study, we extracted DNA from sporocarps consisting of dikaryotic hyphae (n + n).

In Northern Europe, *P. nigrolimitatus* decomposes primarily logs of Norway spruce (Picea abies), but it may also, in rare instance, occur on Scots pine (Pinus sylvestris) (Stokland & Kauserud, 2004). The confounding effects due to different host tree species were avoided by collecting sporocarps exclusively from spruce logs, as the sporocarps on Scots pine might display a different pattern of life-history evolution driven by the host. We collected sporocarps of *P. nigrolimitatus* from 24 spruce-dominated locations in Northern Europe, mainly during autumn 2016 and 2017. Most samples were collected from three regions of Fennoscandia, southeast Norway (NOR), southwest (SW) Finland and northeast (NE) Finland, plus an additional site in central Sweden. Outside Fennoscandia, we sampled one site in northwestern Russia and another in eastern Poland (Table S1, Figure 1). In total, we included 327 sporocarps, collected on independent Norway spruce logs, with a sample size varying of 5–26 per location (Table S1). For dead wood-dependent fungi, the distance between two sampled sporocarps (from different logs) varied according to the availability of dead wood in the forest. As this species typically spreads by sexual basidiospores, no clonal dispersal between logs was expected.

2.2 | Sequencing and annotation of the
reference genome

The dikaryotic strain *P. nigrolimitatus* SigPhenig9 was sequenced and annotated to generate a reference genome for this study. The fruit body was collected from Norway (Sigdal location) and the strain was isolated on malt extract agar (MEA) supplemented with 1 mg/L of benomyl and antibiotics (33 mg/L ampicillin, 8.3 mg/L streptomycin, 3 mg/L tetracycline). For the extractions of DNA and RNA, mycelium was harvested after 20 days of culture in liquid MEA with gentle shaking. The mycelium was filtered and ground in liquid nitrogen. Genomic DNA was extracted in 50 ml Falcon tubes containing 10 ml of lysis buffer (0.2 M Tris pH8.5, 0.25 M NaCl, 25 mM EDTA and 0.5% SDS), followed by a phenol chloroform purification, an RNase treatment and an additional purification using the DNA gTip (Qiagen Genomic-tip 500/G). RNA was extracted with the midi RNA kit (Qiagen) according to the manufacturer's instructions. DNA and RNA were quantified with the Qubit dsDNA BR and Qubit RNA HS Assay kits (Life Technologies), respectively, and visualized by electrophoresis in a 1% agarose gel to check size and quality.

The reference genome was sequenced on Pacific Bioscience (PacBio) platform. We sheared 1 µg of DNA into 10 kb fragments with the g-TUBE(TM) (Covaris) and treated it with DNA damage repair mix, followed by end repair and the ligation of SMRT adapters with the PacBio SMRTbell Template Prep Kit (PacBio). The PacBio sequencing primer was annealed to the SMRTbell template library and the Sequel Binding kit 2.0 was used for sequencing. The library was sequenced on a Pacific Biosystem's Sequel sequencer with the v3 sequencing primer, 1 M v2 SMRT cells, and v2.0 sequencing chemistry, with 6- and 10-hour sequencing run times. Filtered subread data were assembled, with Falcon v1.8.8 (https://github.com/PacificBiosciences/FALCON), and error corrected with Arrow version SMRTLink v5.0.1.9578 (https://github.com/PacificBiosciences/GenomicConsensus).

The transcriptome was sequenced by Illumina RNASeq meth-
ods with PolyA selection. Plate-based RNA sample preparation was

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

FIGURE 1 Maps showing the 24 sampling locations from which dead wood fungus *Phellopilus nigrolimitatus* was collected in spruce-dominated forests in Northern Europe. The colour gradient indicates the mean values for the population genetics statistics measured as expected heterozygosity (a), nucleotide diversity (b) and Tajima's D (c) [Colour figure can be viewed at wileyonlinelibrary.com]
performed on a PerkinElmer Sciclone NGS robotic liquid handling system with the Illumina's TruSeq Stranded mRNA HT sample prep kit and poly-A selection of mRNA according to the Illumina protocol (https://support.illumina.com/sequencing/sequencing_kits/truseq-stranded-mrna.html), with a total starting material of 1 µg of RNA per sample and eight cycles of PCR for library amplification. The prepared library was quantified with the KAPA Biosystem next-generation sequencing library qPCR kit and run on a Roche LightCycler 480 real-time PCR instrument. The quantified library was multiplexed and prepared for sequencing on the Illumina HiSeq sequencing platform with the TruSeq paired-end cluster kit v4 and an Illumina cBot instrument to generate a clustered flow cell for sequencing. The flow cell was sequenced on the Illumina HiSeq 2500 sequencer with HiSeq TruSeq SBS sequencing kits, chemistry v4, according to a 2 × 150 indexed run protocol. Quality-filtered and trimmed RNA-Seq reads were assembled into consensus sequences with Trinity v.2.1.1 (Grabherr et al., 2011). The genome was annotated based on the transcriptome and with the JGI Annotation pipeline (Kuo et al., 2014) to generate a list of 14,792 protein-coding genes and their annotations (Grigoriev et al., 2014).

## 2.3 | DNA extraction and RAD library preparation

Phellopilus nigrolimitatus produces perennial sporocarps, the outermost part of which corresponds to the most recent growth. DNA yield was maximized by cutting small pieces of tissue from the outermost part of the sporocarp and collecting them in 2 ml Eppendorf tubes. The fungal material was ground with a Qiagen Tissue Lyser with 2 mm tungsten beads at 40 Hz for 4 × 45 s. Genomic DNA was extracted with a modified CTAB protocol (Gardes & Bruns, 1993; Murray & Thompson, 1980). We added 600 µl of 2X CTAB and 0.2% β-mercaptoethanol to the sample, which was then incubated at 65°C for 45 min, cooled and mixed with 400 µl chloroform, before centrifuging at 13,000 × g for 10 min. We added an equal volume of cold isopropanol to the supernatant for DNA precipitation overnight at -20°C. The DNA was washed twice with cold 70% ethanol, the extracted gDNA was eluted in 50 µl of elution buffer and cleaned with the E.Z.N.A Soil DNA kit (Omega Biotek) in accordance with the manufacturer's instructions, but with elution in 100 µl EB. We measured the DNA concentration of each sample with the Qubit dsDNA BR Assay kit (Life Technologies) and used 10 mM Tris to standardize the concentration to about 40 ng/µl.

RADseq was performed as described by Maurice et al. (2019). Briefly, the DNA was first digested with the restriction enzyme PstI-HF and ligated to the cutting site with P1 adapters. Twelve uniquely barcoded samples were pooled into each library and sheared by sonication with a 2100 Bioruptor (Diagenode), through four cycles of 30 s ON/30 s OFF. Libraries were concentrated (to 30 µl) with a MinElute Gel Extraction Kit (Qiagen) and subjected to size selection with the Blue Pippin Size Selector (Sage Sciences). Eluted 350-650 bp size-selected libraries were purified with AMPure beads (Agencourt), and ligated to P2 adapters. We amplified

the libraries in a reaction volume of 25 µl and cleaned the PCR products with AMPure beads. The 22 libraries were quantified by qPCR, pooled in equimolar concentrations and subjected to 150 bp paired-end sequencing in two Illumina HiSeq 4000 lanes at the Norwegian Sequencing Centre.

## 2.4 | RADseq data analysis

The sequenced reads were demultiplexed with the *process_radtags* program in *stacks* v.2.10 (Rochette et al., 2019). All reads were truncated to 140 bp and only paired reads were retained for downstream analyses. The demultiplexed reads were mapped to the 141 primary scaffolds of the reference genome of *P. nigrolimitatus* SigPhenig9 (https://mycocosm.jgi.doe.gov/Pheni1) with *BOWTIE2* (Langmead & Salzberg, 2012). The paired-end reads were mapped with a maximum distance of 1000 bp allowed within the pair. The resulting *SAM* files were viewed and ordered in Samtools and converted into *BAM* files. Prior to variant calling, each *BAM* file was assigned a read group and duplicate reads were marked with Picard tools (Broad Institute, 2015). Variants were called with using the *GVCF* workflow in *GATK* v.4.0.0.0 (Poplin et al., 2018). Variants were first called by individual sample with the haplotype caller. The *gvcf* files (containing the individual variant calls) were merged with *CombineGVCFs* and genotyped with *GENOTYPEGVCF*. Variants were excluded with *VariantFiltration* in *GATK* (QD (quality by depth) <2, FS (FisherStrand) >60, MQ (RMSMappingQuality) <40, MQRankSum <-12.5 and ReadPosRankSum <-8). Further filtering of variants and individuals was performed with *VCFTOOLS* v0.1.16 (Danecek et al., 2011), retaining only variants with minGQ (genotype quality) >20. After removing individual samples with high level of missing data (>50%), we then filtered the remaining samples with a maximum missingness of variants threshold at 20%. For the downstream analyses we also removed SNPs with a minimum allele frequency of <0.05, with the --maf function in *VCFTOOLS*.

## 2.5 | Genetic diversity and divergence

For each location sampled, we calculated summary statistics, such as the expected and observed heterozygosity and $F_{IS}$ with *GENODIVE* (Meirmans & Van Tienderen, 2004). Nucleotide diversity indices (theta and $\pi$), and estimates of the departure from neutrality estimate (*Tajima's D*) were calculated with *ANGSD* (Korneliussen et al., 2014). In *ANGSD*, the site frequency spectrum (SFS) for each location was first estimated from the paired-end *BAM* files and the resulting SFS was then used as a prior for the calculation of each statistic per site with the doTheta command. Only SNPs with a coverage above 8 in 50% of the individuals within a location were retained, while SNPs with low quality (minQ < 20 and minMapQ < 30) were excluded from the analysis. We calculated Tajima's D and $\pi$ within 50 kb windows, with a step size of 10 kb, with the thetaStat program in *ANGSD*. Linkage disequilibrium (LD) between loci was estimated with the
geno-r2 command in *VCFTOOLS*. We calculated pairwise LD values at a maximum distance of 50 kb. For the visualization of linkage decay, the mean R² was calculated for pairs of markers with the same physical distance and plotted with R.

Population structure was studied using three different approaches, (i) Bayesian clustering analysis (ADMIXTURE; Alexander et al., 2009), (ii) discriminant analysis of principal components (DAPC; Jombart et al., 2010) and (iii) principal component analysis (PCA). We ran ADMIXTURE with K values from 1-10 with 10 replicates for each of the 10 runs. The results of the replicate runs for each K value were combined, with CLUMPP (Jakobsson & Rosenberg, 2007), and plotted with DISTRUCT (Rosenberg, 2004), using the online CLUMPAK pipeline (Kopelman et al., 2015). We visualized the geographic distribution of the genetic groups identified in ADMIXTURE by plotting the average membership for each sampling location with QGIS v. 2.18.5 (QGIS Development Team, 2018. Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org). We performed a DAPC analysis with the R package ADEGENET (Jombart & Ahmed, 2011), according to the instructions provided in the manual. We used the find.cluster function to identify the number of clusters and included all the principal components (PCs) to capture the maximum variation in the data set. The Bayesian information criterion (BIC) was used to identify the most meaningful number of clusters in the data set. We also applied the dapc function to transform the data by PCA and we performed discriminate analysis (DA) on the retained PCs. The retention of too many PCs can affect the DA, we therefore followed the procedure of Miller et al. (2020). We performed the DAPC in two rounds, in which we first retained 30 PCs and applied the optim.a.score function to identify the optimal number of PCs, which we then used in the second DAPC. For further investigations of population structure, we performed PCA implemented in the R package SNPRELATE, with the dosage mean value used for the impulation of missing data (Zheng et al., 2012).

The degree of genetic differentiation between locations was estimated with the fixation index ($F_{ST}$). $F_{ST}$ values between all pairs of locations were calculated with a custom script in R, using the weighted estimator developed by (Hudson et al., 1992). We assessed isolation by distance (IBD), by calculating the correlation between genetic differentiation (linearized genetic differentiation: $\frac{F_{ST}}{1-F_{ST}}$) and geographic distance in the R package ADEGENET (Jombart & Ahmed, 2011; Rousset, 1997) and estimating the significance of the correlation in a Mantel test (Mantel, 1967), with the mantel.randtest command in the ade4 package (Dray & Dufour, 2007).

## 2.6 | Demographic analysis

We studied the demographic history of *P. nigrolimitatus* in Fennoscandia by comparing three scenarios modeling the split between Finland and Norway, in which samples were collected from multiple locations (unlike Sweden, Russia and Poland in

which samples originated from single sites). Characterization of the split and later gene flow is crucial if we are to understand the colonization of Fennoscandia by *P. nigrolimitatus*. The first scenario modelled a splitting of *P. nigrolimitatus* into two populations (Finland and Norway) at some time in the past, with no contact after the split. In the second scenario, isolation with migration (IM), the two populations were connected by gene flow after the split. In the third scenario, the two populations were isolated for a period of time, before the establishment of a later secondary contact. These three scenarios were investigated by coalescence simulations in FASTSIMCOAL2 (Excoffier et al., 2013; Excoffier & Foll, 2011). The likelihood of each model, and point estimates for demographic parameters, such as the effective population size ($N_e$) of the two populations (logarithms of the two effective population sizes, prior = Unif (2-5)), time since the split (i.e., in generations, prior = Unif (50-300), but note that the upper bound is not a hard constraint in Fastsimcoal2) and the level of bidirectional gene flow (only in models 2 and 3, prior = unif (0.0001-0.1)) were inferred from the simulated site frequency spectrum fitted to the folded observed minor allele frequency spectrum (including invariant sites) with the composite likelihood calculated in FASTSIMCOAL2. For each model, 10 independent FASTSIMCOAL2 runs were performed, each including 1,000,000 coalescence simulations and 60 ECM cycles for likelihood maximization. Confidence intervals for the point estimates were calculated by the parametric bootstrap approach, in accordance with the instructions in the FASTSIMCOAL2 manual. This was achieved by bootstrapping the parameter file with the highest likelihood 1000 times for each model and obtaining the point estimates for the bootstrapped parameter files. Mutation rate is an important parameter determining the shape of the SFS, but has not been precisely estimated for basidiomycetes. Based on both literature and demographic analyses of fungi with different life cycles and hosts, we applied an approximate mutation rate of $1 \times 10^{-7}$.

## 2.7 | Genotype-environment associations

The associations between genotypes and environmental variables were explored by two different approaches, (i) a multivariate analysis for identifying outlier loci in a redundancy discriminate analysis (RDA) and (ii) a univariate logistic regression as implemented in SAMβADA v.0.7 (Stucki et al., 2017). The environmental metadata included the 19 standard bioclimatic variables for Worldclim2 (http://worldclim.org/version2) with a resolution of 30 seconds. We also inferred the level of ammonia NHx, which has been reported to impact the distribution of saprophytic fungal species (Andrew et al., 2018). To reduce correlation among the climatic variables included in the RDA analysis, we excluded all variables with a Spearman's correlation coefficient (|ρ|) >0.7, resulting in a selection of six climatic variables in the RDA analysis: mean annual temperature (Bioclim 01), isothermality (Bioclim 03), maximum temperature of warmest month (Bioclim 05), temperature annual
range (Bioclim 07), mean temperature of driest quarter (Bioclim 09) and annual precipitation (Bioclim 12). The RDA is constrained by the environmental variables while maximizing the amount of genotypic variation explained by these environmental variables. This method allows the identification of outlier loci, which are significantly more strongly associated with the environmental variation constraining the RDA axes than average. The outlier loci are more likely to have been influenced by selection than the other loci, which are more affected by other evolutionary forces, such as genetic drift and gene flow. This method for detecting the genetic loci influenced by selection has been shown to generate very few false-positive signals of selection for populations with a shallow structure (Forester et al., 2018). We performed RDA with the Vegan R package (Oksanen et al., 2019) and identified outlier loci on all significant RDA axes. All loci with a loading exceeding three standard deviations from the mean were considered significantly associated with environment (Forester et al., 2018).

In SAMβADA, we calculated the univariate association of each SNP genotype, coded as present or absent, with the set of 20 environmental variables. Despite collinearity, we retained all variables in the SAMβADA analyses, which are univariate, implying that the analysis of one variable does not affect the result for another variable. The resulting *p*-values were adjusted for false discovery rate, by calculating *q* values in the R package *QVALUE* (Storey et al., 2019), and the genotype-environmental associations with a *q*-value below 0.05 were retained.

To visualize the potential effect of the environmental associations on the allele frequencies and the pattern of the associations, we calculated the Spearman's correlation coefficients between the allele frequencies per population from each SNP significantly association with the environmental variable in either SAMβADA or the RDA, using the rcorr function in the Hmisc R package (Harrell, 2021). We calculated the Euclidean distance between SNPs based on their correlations with environmental variables, using the dist function, and the distances obtained were used to cluster the SNPs with the hclust function in R. We visualized these correlations with the heatmap.2 function in the GPLOTS package (Warnes et al., 2020).

## 2.8 | Annotation of loci under selection

The resulting SNPs significantly associated with environmental variables were annotated with the JGI gene repertoire and genome assembly for *P. nigrolimitatus*. We identified the genes overlapping the SNPs or the first genes downstream from the SNPs that were significantly associated with climatic variables using the closest command in BEDTOOLS (Quinlan & Hall, 2010). We further investigated whether there was any significant enrichment in protein domains (IPR) among the loci associated with climatic variation. The significance of the enrichment was determined by performing a Fisher's exact test in R and the *p*-values were adjusted for multiple tests with the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995), with the p-adjust command in R.

# 3 | RESULTS

## 3.1 | The reference genome

We obtained a reference genome for subsequent analyses by sequencing the dikaryotic isolate SigPhenig9 of *Phellopilus nigrolimitatus* on the Pacific Biosciences sequencing platform. The genome was assembled into 696 scaffolds with a total length of 56.76 Mb. When two haplotypes of a dikaryotic genome are divergent, they produce a primary and a secondary scaffold that are homologous. Exclusion of the homologous secondary scaffolds resulted in 141 primary scaffolds with a total length of 39.68 Mb, with scaffold N50 of 17 and L50 of 640 Kb (see Table S2 for the genome statistics calculated for primary scaffolds only). The remaining 555 scaffolds were highly similar to larger scaffolds, and therefore probably represented allelic variants. From the primary scaffolds, we predicted the presence of 14,792 genes, including ~2000 members of multigene families, the largest of which included genes encoding protein kinases (283 copies), cytochrome P450 (144 copies), and major facilitator superfamily transporters (133 copies). The genome also encodes for 72 haloacid dehalogenase-like hydrolases, one of the largest numbers of such enzymes in wood-decay fungi.

## 3.2 | Genetic diversity

RAD sequencing of samples from the 24 locations provided a uniform coverage of the *P. nigrolimitatus* primary scaffolds, spanning a mean of 25% of the genome (Figure S1 Distribution of SNPs across all primary scaffolds), with a mean coverage of 23.75 reads per individual. After the exclusion of individuals and genomic sites for which read quality was insufficient, we retained 93,938 biallelic SNPs from 279 individuals for downstream analyses.

Mean nucleotide diversity ($\pi$) was $8.52 \times 10^{-3}$ per kb and ranged from $7.98 \times 10^{-3}$ and $9.32 \times 10^{-3}$ per kb at the various locations (Table S1, Figure 1). Mean expected heterozygosity ($H_e$) was 0.103 and ranged from 0.097 to 0.108. Southwestern (SW) Finland harbored the lowest diversity, for both $\pi$ (mean = $8.34 \times 10^{-3}$) and $H_e$ (mean = 0.101), but genetic diversity values were generally higher in northeastern (NE) Finland (mean $\pi$ = $8.65 \times 10^{-3}$ and mean $H_e$ = 0.105) and Norway (mean $\pi$ = $8.70 \times 10^{-3}$ and mean $H_e$ = 0.104, see Figure 1). We assessed deviation from the neutral model of molecular evolution, by calculating Tajima's D, which was negative in all locations, and ranged from -1.317 to -0.047. In general, the most negative Tajima's D values were found in the Norway (for all but the most western site, Mørkvasslia), Sweden and Poland, whereas values in SW Finland were close to zero and those in NE Finland were intermediate (Table S1, Figure 1).

Linkage disequilibrium (LD) decayed rapidly, whether estimated from total samples or assessed regionally, implying high rates of recombination (Figure S2 Linkage disequilibrium decay). In all calculations, the mean LD for all samples was below 0.1 at a distance of 100 bp. When LD was calculated within regions (Norway, SW
Finland, NE Finland) a similar trend was observed, with rapid LD decay in all regions. LD values were slightly higher in Norway and NE Finland than in SW Finland, which is surprising as LD is generally higher in smaller populations, such as those in SW Finland. However, the general linkage disequilibrium pattern suggests that processes such as genetic drift and inbreeding, which could increase the genome-wide LD, play only a minor role in structuring the genetic variation in *P. nigrolimitatus*.

## 3.3 | Genetic divergence

An analysis of the population structure of *P. nigrolimitatus* with ADMIXTURE showed that the lowest cross-validation error was obtained for *K* = 1, suggesting a very shallow population structure (Figure S3). For *K* = 2, the Norwegian and the Finnish samples clustered into two separate groups (Figure 2). A similar pattern was obtained for *K* = 3, with an additional group consisting of three locations from SW Finland. These locations — Luukki, Herukkapuro and Vesijako — were among those with the lowest genetic diversity and the highest Tajima's D. For *K* = 4, the Finnish locations were split into a SW Finnish and NE Finnish group, Sweden was intermediate between Norway and Finland, and the Russian and Polish locations were more similar to NE Finland and SW Finland, respectively. The structure for higher *K* values is presented in Figure S4. The DAPC analysis

also identified *K* = 1 as the most likely structure (Figure S5). When analysing the structure at *K* = 2, the second group consisted of only four individuals from Vasijako. For larger cluster numbers, the structure obtained with the DAPC analysis was similar to that suggested by principal component analysis (PCA, Figure S6). The PCA did not distinctly separate Norwegian from Finnish locations, and the first two axes of the PCA explained very low proportions of the variance: 0.90% for axis 1 and 0.85% for axis 2. PCs 1 and 2 separated the three SW Finnish locations forming a separate group in the ADMIXTURE analysis (*K* = 3, Figure 2b). The exclusion of these locations from the PCA resulted in a slight differentiation between the Norwegian and Finnish locations, with the Swedish location in an intermediate position (Figure S6).

The PCA results were supported by the relatively low pairwise *F*<sub>ST</sub> estimates, ranging between 0 and 0.05 (Figure 3a), except for the locations in SW Finland that also deviated in the PCA and ADMIXTURE analyses. There was no significant isolation by distance in the analysis of all locations (Figure 3b), probably due to the SW Finnish locations displaying high levels of genetic differentiation from geographically close locations. The opposite trend was observed for Poland, with very low genetic differentiation despite the relatively large geographic distance from other locations. However, the exclusion of all locations in SW Finland and Poland from the data set, revealed significant isolation by distance (*r* = 0.27, *p*-value = .013), probably reflecting the natural dispersal of *P. nigrolimitatus* after the last period of glaciation.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Admixture analysis of *Phellopilus nigrolimitatus*. The mean individual membership ancestry coefficient per population is plotted on the map for *K* = 2 (a), *K* = 3 (b) and, *K* = 4 (c). Individual membership is plotted for *K* = 2–4 (d) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Heat map of pairwise $F_{ST}/(1-F_{ST})$ values between populations of *Phellopilus nigrolimitatus*. Dark colour with high $F_{ST}/(1-F_{ST})$ values indicates a higher level of genetic differentiation. (a). Isolation by distance (IBD) in *Phellopilus nigrolimitatus*. Grey circles and line indicate the isolation by distance for all 24 populations included, whereas the red circles and line indicate isolation by distance excluding the populations from Poland and southwestern Finland (b) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<table><caption>TABLE 1 Demographic inference of the fungus *P. nigrolimitatus* in Fennoscandia. Point estimates generated with Fastsimcoal2 from three demographic models (1) population split between Norway and Finland, (2) isolation by migration with gene flow, and (3) secondary contact after isolation</caption><thead><tr><th>Model</th><th>N<sub>e</sub> (Fin)</th><th>N<sub>e</sub> (Nor)</th><th>T<sub>div</sub></th><th>T<sub>div2</sub></th><th>m<sub>N-F</sub></th><th>m<sub>F-N</sub></th><th>AIC</th></tr></thead><tbody><tr><td>Divergence without migration</td><td>8055 (7989.2-8109.6)</td><td>8053 (7504.6-8273.6)</td><td>107.3 (104.1-111.1)</td><td>NA<sup>a</sup></td><td>NA</td><td>NA</td><td>259424.9</td></tr><tr><td>Divergence with migration</td><td>7158 (7152.6-7312)</td><td>8442 (8150.2-9544.2)</td><td>296.0 (258-296)</td><td>NA</td><td>0.0015 (0.0008-0.0018)</td><td>0.0036 (0.0031-0.0038)</td><td>258983.8</td></tr><tr><td>Secondary contact</td><td>7246 (6473.5-7256.0)</td><td>8760 (7654.4-10938.0)</td><td>296.4 (240.9-312.6)</td><td>6.67 (1.61-55.05)</td><td>0.0224 (0.002-0.102)</td><td>0.0589 (0.008-0.254)</td><td>258956.3</td></tr></tbody></table>

Note: In these models, $N_e$ is the effective population size in Norway (Nor) and Finland (Fin), $T_{div}$ is the divergence time (number of generations between Norway and Finland), $m_{N-F}$ is the rate of migration of individuals from Norway to Finland and $m_{F-N}$ is the rate of migration from Finland to Norway. Value ranges in parentheses indicate confidence intervals.

<sup>a</sup>NA indicates parameters that were not estimated in some scenarios.

## 3.4 | Historical demography

The two demographic scenarios including gene flow between the Norwegian and the Finnish populations (model 2 and 3) had the highest likelihood (Table 1, Figure S7 Demographic models). The key difference between these scenarios is continuous gene flow after the splitting of the populations in scenario 2, whereas scenario 3 includes a separation of the populations for a period of time before the occurrence of secondary contact leading to gene flow. The difference in the likelihoods of these two scenarios was small, and their demographic estimates were similar. Both scenarios suggest a separation about 296 generations ago (Table 1), together with asymmetric

migration, with greater gene flow from Finland into Norway than in the opposite direction. However, gene flow rates are much higher in the scenario involving secondary contact than in the IM scenario, in which gene flow is assumed to have occurred only in the last six to seven generations.

## 3.5 | Genotype - environment associations

The RDA analysis revealed a stronger differentiation of the locations along the axes (Figure 4a) constrained by climatic variables than along the PCA axes (Figure S6). Though only, a very small proportion
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Redundancy analyses of the SNP variation explained by environmental variables. The points indicate the 279 individuals of *P. nigrolimitatus* (coloured by geographic region) inferred in the RDA analyses and projected along the first two axes RDA1 and RDA2, which explain 0.7% and 0.5% of the variation respectively (a). Projection of 27 607 SNPs and six explanatory environmental variables (|r| > 0.7) in the RDA. Of the 836 significant SNPs, 249 (blue) were specific to the RDA, 437 (green) were identified in the SAMβADA analysis and 150 (yellow) were common to SAMβADA and RDA (b) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

of the variation can be explained by the RDA axes, there are nev-
ertheless, associations between genotype and environment in *P.
nigrolimitatus*. The first five RDA axes explained 3.0% of the total
genetic variation and we identified 399 SNPs as outliers (Figure 4b
and Figure S8).

We further identified 587 SNPs significantly associated with at least one environmental variable, with *SAMβADA*. Given the correlations between several environmental variables, many SNPs were significantly associated with more than one variable (Figure S9 shows the number of significant tests per environmental variable). Figure 5 shows the Manhattan plots of the *p*-values for the six environmental variables selected in the RDA along the scaffolds. A clear pattern is visible for each variable, with signals of selection detected in different genomic regions. In total, 150 of the 587 SNPs significantly associated with environmental variables in *SAMβADA* overlapped with the findings of the RDA (Figure 4b). In total, 836 SNPs were associated with climatic variables in one of

the two analyses. Furthermore, the SNPs identified as significantly
associated with climatic variables in SAMβADA were also among
those with the strongest associations in the RDA (along RDA axes
1 and 2 in particular; Figure 4). Many of the SNPs associated with
climatic variables were located within exons (568), 120 were in in-
trons within genes, and the remaining 148 SNPs were located in
intergenic regions.

There are several explanations for the finding of significant asso-
ciations with genomic regions outside genes. For example, significant
SNPs in intergenic regions may be linked to functional variation in
genes not targeted by the RAD sequencing approach, or may be lo-
cated in regions important for the regulation of one or several genes.
The 568 intragenic SNPs were present in a total of 340 genes, 261 of
which had functional annotations (the remaining 79 were considered
to correspond to hypothetical proteins). The strongest association in
the SAMβADA analysis was that for two SNPs (position 342,225 and
342,533) in scaffold 34. These SNPs were tightly linked to a 4 bp
<figure><img src="image_6.png" /></figure>

FIGURE 5 Manhattan plot of the *p*-values (-log₁₀) over the genome of *P. nigrolimitatus* for the loci associated with six selected environmental variables from Worldclim2: (a) annual mean temperature, (b) isothermality, (c) maximum temperature of warmest month, (d) temperature annual range, (e) mean temperature of driest quarter, and (f) annual precipitation. Each point represents a genotype-environment association and the points above the horizontal line (red) are significant after adjustment for false discovery rate (q-value < 0.05) [Colour figure can be viewed at wileyonlinelibrary.com]

indel (excluded from the analysis) and one of the SNPs (342,533) was located in a coding region. This genomic region (scaffold 34) is very similar to another region in scaffold 10, possibly resulting from an earlier duplication event. No annotation was available for this region. We therefore performed a Blastp search for the translated sequence in the coding region (34,205–342,598) to identify genomic regions encoding similar amino-acid sequences. This regulator has not yet been found in any eukaryote. The only similar sequence identified to date was that of an Lrp/AsnC family transcriptional regulator that is conserved in bacteria and interestingly in *Collimonas fungivorans*, a mycophagous bacterium with ability to grow and feed on fungi (Mela et al., 2011).

To identify the overall patterns in genotype-environmental asso-
ciations, we extracted all 836 SNPs significantly associated with any
environmental variable and performed a hierarchical cluster analy-
sis of the Spearman's rank correlation coefficients between allele
frequencies and environmental variables. This analysis identified
two principal groups of SNPs, with related correlations with envi-
ronmental variables (Figure 6a). The first group contained 304 SNPs
mostly correlated with temperature variables (both mean annual

temperature and temperature in the coldest periods) and season-
ality. The second group of 532 SNPs was more strongly associated
with temperature in the wettest and warmest quarters. The correla-
tion between SNPs and environmental variables followed a very
similar pattern as the correlation between the different environmen-
tal variables (Figure 6b), suggesting that the correlation of several
environmental variables with the same SNP probably reflects the
autocorrelation of environmental variables rather than environmen-
tal pleiotropy.

We performed a functional enrichment analysis for annota-
tions of interprotein domains in genes containing SNPs associ-
ated with environmental variables. We first investigated the set
of 150 SNPs that was significant in both RDA and SAMβADA and
identified several functional domains displaying significant en-
richment (Table 2). The functional groups for which enrichment
was strongest were the F-box domains, small GTPases and heli-
cases. We performed a similar enrichment analysis on all SNPs
associated with environmental variation in either RDA or SAMβADA.
We performed separate analyses for the two different groups
identified in the cluster analysis (see Figure 6a) and found some
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Associations between environmental variables and SNPs. Hierarchical clustering of associations between SNPs (columns) and environmental variables (rows) assessed by calculating Spearman's ρ in absolute values suggest two main clusters (groups 1 and 2) (a). Correlation between environmental variables assessed by calculating Spearman's ρ according to the hierarchical clustering of environments, revealed a pattern similar to that for SNP-environment associations (b) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<table><caption>TABLE 2 Interprotein domains (Interprotein) significantly enriched in the SNPs that associated with climate in both RDA and SAMβADA</caption><thead><tr><th>Interprotein</th><th>Function</th><th>p-value<sup>a</sup></th><th>Freq<sup>b</sup></th><th>Freq_all<sup>c</sup></th></tr></thead><tbody><tr><td>IPR001810</td><td>F-box domain</td><td>***</td><td>7</td><td>67</td></tr><tr><td>IPR005225</td><td>Small GTP-binding protein domain</td><td>***</td><td>3</td><td>10</td></tr><tr><td>IPR001806</td><td>Small GTPases</td><td>**</td><td>2</td><td>4</td></tr><tr><td>IPR003578</td><td>Small GTPase Rho</td><td>**</td><td>2</td><td>4</td></tr><tr><td>IPR003579</td><td>Ras small GTPase</td><td>**</td><td>2</td><td>4</td></tr><tr><td>IPR020849</td><td>Small GTPase superfamily</td><td>**</td><td>2</td><td>4</td></tr><tr><td>IPR000008</td><td>C2 domain</td><td>**</td><td>2</td><td>7</td></tr><tr><td>IPR000629</td><td>ATP-dependent RNA helicase DEAD-box, conserved site</td><td>**</td><td>2</td><td>7</td></tr><tr><td>IPR009000</td><td>Translation protein, beta-barrel domain superfamily</td><td>*</td><td>2</td><td>11</td></tr><tr><td>IPR014014</td><td>RNA helicase, DEAD-box type, Q motif</td><td>*</td><td>2</td><td>12</td></tr><tr><td>IPR027417</td><td>P-loop containing nucleoside triphosphate hydrolase</td><td>*</td><td>8</td><td>182</td></tr><tr><td>IPR001650</td><td>Helicase, C-terminal</td><td>*</td><td>3</td><td>41</td></tr><tr><td>IPR014001</td><td>Helicase superfamily 1/2, ATP-binding domain</td><td>*</td><td>3</td><td>41</td></tr></tbody></table>

differences in enrichment in functional domains between the two groups of SNPs with different environmental associations (Table S3). The first group was mostly enrichment in GTPases and Ankryn

repeats, whereas the second was enriched in GTPases, helicases and heavy metal associated domains, and also cellulose/chitin-binding proteins.

<footnote><sup>a</sup>***p-value <.001, **p-value <.01 *p-value <.05. p-values in bold are significant after correction for multiple testing.</footnote>

<footnote><sup>b</sup>Frequency (Freq.) of genes from an interprotein domain in which a significant SNP was found.</footnote>

<footnote><sup>c</sup>Freq_all represent the frequencies of each interprotein domain in the genome.</footnote>
4 | DISCUSSION

**4.1 | Genetic diversity and population structure**

Several lines of evidence suggest that the observed population ge-
netic structure of *P. nigrolimitatus* results from a combination of gene
flow, genetic drift and selection. *Phellopilus nigrolimitatus* probably
followed the westward expansion of its main host, Norway spruce,
after the last glaciation period. On a regional scale, *P. nigrolimitatus* is
characterized by a general trend towards weak genetic structuring,
consistent with earlier population genetic studies of other wood-
decay fungi in Fennoscandia (Franzén et al., 2007; Högberg et al.,
1999; Kauserud & Schumacher, 2002). Higher divergence of some
populations in SW Finland, along with a more recent decline in popu-
lation size, may be the consequence of fragmentation and loss of
old-growth forests at the landscape-level.

A separation of the Norwegian populations from those of NE
Finland and Russia, is in line with the genetic structure of Norway
spruce (see Figure 2; Tollefsrud et al., 2009). Yet, genetic differ-
entiation is less pronounced in *P. nigrolimitatus* than its host, sug-
gesting that the fungal populations are connected by higher gene
flow than its host populations. The low FST values between south-
ern Norway and several Finnish locations point to a recent diver-
gence of *P. nigrolimitatus* in Fennoscandia, which is also supported
by coalescence simulations. The recent divergence identified here,
along with gene flow at the regional scale, probably resulted in the
weak genetic structure as revealed by ADMIXTURE, DAPC and PCA. We
estimated the divergence time between Norway and Finland to be
296 generations. The generation time of wood-inhabiting fungi is
largely unknown. *P. nigrolimitatus* has not been found on living trees
nor is known to be able to infect them. However, it has been de-
tected in sawdust drilled from downed logs at relatively early stages
of decay (Ovaskainen et al., 2013), and the species mainly fruits on
old, well-decayed logs (Stokland & Kauserud, 2004). Spruce logs may
take more than 40 years to reach such a level of decay (Mäkinen
et al., 2006), suggesting that the generation time of *P. nigrolimita-
tus* may be between 20–40 years, allowing establishment in already
decayed logs and overlapping generations. Converting the number
of generations into years, with a generation time of 20–40 years,
we speculate that the *P. nigrolimitatus* populations separated about
6000–12,000 years ago (without overlapping generations). This es-
timate corresponds to a separation after the last period of glaciation
and is in line with a split into different genetic clusters identified in
Norway spruce (Tollefsrud et al., 2009). The negative Tajima's D val-
ues obtained for all locations are also indicative of range expansion,
most probably after the last glaciation period. The most negative
values were found at Norwegian locations, at which Norway spruce
has been a dominant tree species for only the last 2000–3000 years
(Tollefsrud et al., 2008).

Land use, including forestry, agriculture and urbanization,
has been more intensive in SW Finland than in the NE Finland
(Kalliola, 1966; Lilja & Kuuluvainen, 2005; Rouvinen et al., 2002).
Southwestern Finland has experienced permanent settlement since

the Stone Age, whereas there have been substantial permanent set-
tlements in NE Finland for only about 300 years (Koskinen et al.,
2007). Settlement history, with human population density and forest
use increasing over time, has resulted in lower proportions of old-
growth forests in SW Finland (Punttila & Ihalainen, 2006), account-
ing for the lower genetic diversity and higher genetic differentiation
at locations in this region. A similar pattern was reported for another
wood-decay fungus, *Phlebia centrifuga*, with displays higher genetic
differentiation in patchy habitats in southern Sweden and lower lev-
els of genetic differentiation between well-connected populations
in Fennoscandia (Franzén et al., 2007). In a meta-analysis, Schlaepfer
et al. (2018) found that nonwoody plants and animals within forest
habitats experienced negative effects of anthropogenic habitat frag-
mentation. Most importantly, the time since habitat fragmentation
was found to affect genetic diversity, with significant negative ef-
fects on habitat fragments that were more than 50 years old show-
ing. Given the long generation time of *P. nigrolimitatus*, one plausible
explanation for the lack of genetic effects of habitat loss in Norway,
NE Finland and NW Russia, is the relatively recent and less extensive
loss and fragmentation of old-growth forests in these areas, result-
ing in the maintenance of habitat connectivity and facilitating gene
flow. A similar situation may exist in Poland, where the analysed
population was genetically similar to those at other locations despite
separation by large geographic distances.

**4.2 | Genotype - environment associations and local adaptation**

Surprisingly little is known about how wood-inhabiting fungi adapt
to climate, despite their ecological and economic importance in the
forest ecosystem. Previous studies of genotype-environment asso-
ciations in fungi have focused on other groups of fungi such as yeasts
(Peter et al., 2018), plant pathogens (e.g., Bartoli & Roux, 2017) and
ectomycorrhizal fungi (Branco et al., 2015, 2017). We hypothesized
that *P. nigrolimitatus* had adapted to different environments dur-
ing and after the westward expansion of its range, but the genetic
variants significantly associated with environmental variables was
surprisingly high. (Branco et al., 2015) identified signals of positive
selection between recently diverged mountain and coastal popula-
tions of *Suillus brevipes*, of which were partly explained by the ad-
aptation to salt tolerance. As for other organisms, there is evidence
that relatively strong local adaptation may evolve within only a few
generations, even in species with high levels of gene flow, like many
forest trees (Eckert et al., 2015; Sork et al., 2016). Long-lived wood-
decay fungi are comparable to trees in several ways. They are sessile
and the mycelium persist over relatively long periods of time. Many
wood-inhabiting fungi also have large populations, rendering selec-
tion more important than genetic drift as an evolutionary force.

We observed a strong association between several climate vari-
ables and two SNPs on scaffold 34. A short coding sequence with
similarity to the Lrp/AsnC family transcriptional regulator in bacteria
was identified in this genomic region. This regulator has not been
identified in any eukaryotes before and, given the low sequence similarity to the Lrp/AsnC family transcriptional regulator, it seems likely that this genomic region has a different function in *Phellopilus nigrolimitatus*. The scaffold 34 genomic region, which is strongly associated with climatic variables, may also be a regulatory region, such regions having often been shown to have strong signatures of selection in other eukaryotes (Jones et al., 2012).

RAD sequencing is a reduced representation of the genome. It cannot, therefore, take into account functional loci strongly associated with environmental variables that were not sequenced. It is thus difficult to identify the actual gene under selection. However, RAD sequencing was successfully used to identify adaptive loci involved in the transition from marine to freshwater life style in threespine stickleback (Hohenlohe et al., 2010). Smith et al. (2019) studied the adaptation of leafcutter ant *Atta texana* and its fungal symbionts using double digested RAD sequencing. In line with other studies of adaptation in fungi, we identified genes involved in membrane transport and secretion (Balasundaram et al., 2018; Branco et al., 2017; Ellison et al., 2011) and helicases known to be associated with environmental factors in the genera *Suillus* and *Neurospora* (Branco et al., 2017; Ellison et al., 2011). Group 1 genes (Figure 6) were associated with environmental stress and were highly correlated with seasonality and temperature and precipitation during the coldest and driest quarter and month. Seasonal differences may require a tight regulation of hyphal growth and secretion. By contrast, group 2 genes, were associated with temperature and precipitation during the warmest month and quarter, potentially reflecting the seasonal timing of events in the fungal life cycle, such as optimum vegetative growth and fruiting. Overall, these examples highlight the importance of diverse regulatory processes in the adaptation of fungi to climate. However, the incomplete annotation of many genes, as is often the case for non-model fungi, may have resulted in gene sets with important functions being overlooked.

5 | CONCLUSION AND PERSPECTIVES

We sequenced and annotated the genome of a wood-inhabiting fungus, *P. nigrolimitatus* and assessed the genome-wide genetic diversity prevailing in 24 populations across the boreal forest of Northern Europe, to gain insights into the genetic basis of fungal adaptation to climate. We found lower levels of genetic diversity in small isolated populations in SW Finland than in the populations of the larger connected landscapes in Northern Fennoscandia, consistent with forest history. However, beyond the landscape of SW Finland, *P. nigrolimitatus* seems to be resilient to changes in genetic diversity, despite being rare in most of these forests. Its life cycle, including long generation time and a high potential for gene flow may contribute to the maintenance of high levels of genetic diversity in *P. nigrolimitatus*. The relatively long generation time of *P. nigrolimitatus* may also buffer the long-term effect of forest fragmentation and consequently slowing the adaptive response of the species to changes in the climate. The diversity at adaptive loci may be of more important information than

neutral genetic diversity in taking conservation decision-making for rare species (Grummer et al., 2019; Lanes et al., 2018). Overall, these results suggest that local adaptation has resulted in the evolution of spatial differences in certain genomic regions.

To conclude, the population genomic analysis of *P. nigrolimitatus* provided new insight into the extent of local adaption in wood decay fungi, which can be compared with the observed adaptation of their well-studied tree host species. The consequences of the combined effect of future climate change and habitat fragmentation of old-growth forests on the loss of biodiversity will be determined by the ability of many forest species to migrate or adapt to new conditions. Knowledge about the adaptive genetic variation, along with climate models, can be included in forest management strategies to predict whether populations have the standing genetic variation that are needed to adapt with future changes in climate (Rochat et al., 2021). However, to fully understand the adaptive potential, we need to experimentally investigate the link between genetic variation and fitness under future environmental conditions. Compared to many other organisms, wood-inhabiting fungi, with the example of *P. nigrolimitatus*, can be isolated and cultured to perform long-term experimental studies that include measures of genotype, phenotype and environmental variation for investigating this correlation.

ACKNOWLEDGEMENTS

We thank Osmo Heikkala, Teppo Helo, Katja Juutilainen, Dariusz Karasinski, Juha Karvonen, Salla Mäkelä, Rose-Marie Oksanen, Jorma Pennanen, Tord Snäll, Irina Sorokina, Viacheslav Spirin and Sten Svantesson for assistance in sample collection. We thank Carrie J. Andrew for providing NHx data. The sequencing, assembly, and annotation of the reference genome were performed by the U.S. Department of Energy Joint Genome Institute, a DOE Office of Science User Facility, supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank Anna Lipzen for the transcriptome assembly, Chris Daum's group for DNA and RNA sequencing of the reference genome and Kerrie Barry for coordinating sequencing. The RAD libraries were sequenced by the Norwegian Sequencing Center hosted by the University of Oslo and supported by the "Functional Genomics" and "Infrastructure" programs of the Research Council of Norway and the south-eastern Regional Health Authorities. SM thanks the group of Guido Barbugiani, University of Ferrara, for hosting her research stay funded by a Kristine Bonnevies grant from the University of Oslo. This work was funded by the Research Council of Norway (grant number 254746). The authors thank Tatiana Giraud and four anonymous reviewers for providing constructive comments on earlier drafts of the manuscript.

AUTHOR CONTRIBUTIONS

Jenni Nordén, Inger Skrede, Håvard Kauserud and Sundy Maurice applied for funding. Jenni Nordén, Håvard Kauserud and Sundy Maurice organized sampling and designed the study. Otto Miettinen contributed with samples. Jørn Henrik Sønstebø, Jenni Nordén and Sundy Maurice conceived the study. Sundy Maurice isolated
reference strain and performed DNA/RNA extractions for WGS. Igor V. Grigoriev and Francis Martin coordinated the sequencing of the reference genome; Jasmyn Pangilinan and Sajeet Haridas assembled and annotated the reference genome. Jørn Henrik Sønstebø led the bioinformatics and statistical analyses and wrote a first draft. Emiliano Trucchi contributed in the demographic analyses. Jørn Henrik Sønstebø and Sundy Maurice performed laboratory work, prepared figures, and wrote the manuscript with contributions from Emiliano Trucchi, Jenni Nordén, Inger Skrede, Håvard Kauserud, Otto Miettinen and Igor V. Grigoriev.

## OPEN RESEARCH BADGES

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

This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.7d7wm37wr.

## DATA AVAILABILITY STATEMENT

The genome sequence of *Phellopilus nigrolimitatus* and its annotations are available from https://mycocosm.jgi.doe.gov/Pheni1. This Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession JAJOHG000000000. The version described in this paper is version JAJOHG010000000. The Fastq files and BAM files generated from the mapping to the P. nigrolimitatus genome are available from the European Nucleotide Archive (project PRJEB48923). The filtered VCF file, environmental metadata and customs scripts used for the population genetic analyses are available from Dryad (doi https://doi.org/10.5061/dryad.7d7wm37wr). Sporocarps and gDNA are available upon request.

## ORCID

Jørn Henrik Sønstebø https://orcid.org/0000-0002-3204-1283

Emiliano Trucchi https://orcid.org/0000-0002-1270-5273

Jenni Nordén https://orcid.org/0000-0001-8894-5815

Inger Skrede https://orcid.org/0000-0002-1113-7403

Otto Miettinen https://orcid.org/0000-0001-7502-710X

Sajeet Haridas https://orcid.org/0000-0002-0229-0975

Jasmyn Pangilinan https://orcid.org/0000-0001-7966-3496

Igor V. Grigoriev https://orcid.org/0000-0002-3136-8903

Francis Martin https://orcid.org/0000-0002-4737-3715

Håvard Kauserud https://orcid.org/0000-0003-2780-6090

Sundy Maurice https://orcid.org/0000-0002-5376-0981

## REFERENCES

Abrego, N., Christensen, M., Bässler, C., Ainsworth, A. M., & Heilmann-Clausen, J. (2017). Understanding the distribution of wood-inhabiting fungi in European beech reserves from species-specific habitat models. *Fungal Ecology, 27*, 168-174. https://doi.org/10.1016/j.funeco.2016.07.006

Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. *Genome Research, 19*(9), 1655-1664. https://doi.org/10.1101/gr.094052.109

Andrew, C., Halvorsen, R., Heegaard, E., Kuyper, T. W., Heilmann-Clausen, J., Krisai-Greilhuber, I., Bässler, C., Egli, S., Gange, A. C., Høiland, K., Kirk, P. M., Senn-Irlet, B., Boddy, L., Büntgen, U., & Kauserud, H. (2018). Continental-scale macrofungal assemblage patterns correlate with climate, soil carbon and nitrogen deposition. *Journal of Biogeography*, 45(8), 1942-1953. https://doi.org/10.1111/jbi.13374

Balasundaram, S. V., Hess, J., Durling, M. B., Moody, S. C., Thorbek, L., Progida, C., LaButti, K., Aerts, A., Barry, K., Grigoriev, I. V., Boddy, L., Högberg, N., Kauserud, H., Eastwood, D. C., & Skrede, I. (2018). The fungus that came in from the cold: Dry rot's pre-adapted ability to invade buildings. *ISME Journal*, 12(3), 791-801. https://doi.org/10.1038/s41396-017-0006-8

Bartoli, C., & Roux, F. (2017). Genome-wide association studies in plant pathosystems: Toward an ecological genomics approach. *Frontiers in Plant Science*, 8, 763. https://doi.org/10.3389/fpls.2017.00763

Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. *Journal of the Royal Statistical Society: Series B (Methodological)*, 57(1), 289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

Boddy, L. (2016). Chapter 4 - genetics - variation, sexuality, and evolution. In S. C. Watkinson, L. Boddy, & N. P. Money (Eds.), *The fungi* (3rd ed., pp. 99-139). Academic Press.

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.

Branco, S., Gladieux, P., Ellison, C. E., Kuo, A., LaButti, K., Lipzen, A., Grigoriev, I. V., Liao, H.-L., Vilgalys, R., Peay, K. G., Taylor, J. W., & Bruns, T. D. (2015). Genetic isolation between two recently diverged populations of a symbiotic fungus. *Molecular Ecology*, 24(11), 2747-2758. https://doi.org/10.1111/mec.13132

Broad Institute (2015). *Picard Tools - By Broad Institute*. Available from, http://broadinstitute.github.io/picard

Chen, J., Li, L., Milesi, P., Jansson, G., Berlin, M., Karlsson, B. O., Aleksic, J., Vendramin, G. G., & Lascoux, M. (2019). Genomic data provide new insights on the demographic history and the extent of recent material transfers in Norway spruce. *Evolutionary Applications*, 12(8), 1539-1551. https://doi.org/10.1111/eva.12801

Dal Grande, F., Sharma, R., Meiser, A., Rolshausen, G., Büdel, B., Mishra, B., Thines, M., Otte, J., Pfenninger, M., & Schmitt, I. (2017). Adaptive differentiation coincides with local bioclimatic conditions along an elevational cline in populations of a lichen-forming fungus. *BMC Evolutionary Biology*, 17(1), 93. https://doi.org/10.1186/s12862-017-0929-8

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S. T., McVean, G., Durbin, R. & 1000 Genomes Project Analysis Group (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156-2158. https://doi.org/10.1093/bioinformatics/btr330

Dray, S., & Dufour, A.-B. (2007). The ade4 Package: Implementing the duality diagram for ecologists. *Journal of Statistical Software*, 22(4), 1-20. https://doi.org/10.18637/jss.v022.i04

Eckert, A. J., Maloney, P. E., Vogler, D. R., Jensen, C. E., Mix, A. D., & Neale, D. B. (2015). Local adaptation at fine spatial scales: An example from sugar pine (*Pinus lambertiana*, Pinaceae). *Tree Genetics & Genomes*, 11(3), 42. https://doi.org/10.1007/s11295-015-0863-0

Ellison, C. E., Hall, C., Kowbel, D., Welch, J., Brem, R. B., Glass, N. L., & Taylor, J. W. (2011). Population genomics and local adaptation in wild isolates of a model microbial eukaryote. *Proceedings of the National Academy of Sciences*, 108, 2831-2836. https://doi.org/10.1073/pnas.1014971108

Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V. C., & Foll, M. (2013). Robust demographic inference from genomic and SNP
data. *Plos Genetics*, 9(10), e1003905. https://doi.org/10.1371/journal.pgen.1003905

Excoffier, L., & Foll, M. (2011). fastsimcoal: A continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios. *Bioinformatics*, 27(9), 1332-1334. https://doi.org/10.1093/bioinformatics/btr124

Forester, B. R., Lasky, J. R., Wagner, H. H., & Urban, D. L. (2018). Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. *Molecular Ecology*, 27(9), 2215-2233. https://doi.org/10.1111/mec.14584

Franzén, I., Vasaitis, R., Penttilä, R., & Stenlid, J. (2007). Population genetics of the wood-decay fungus *Phlebia centrifuga* P. Karst. in fragmented and continuous habitats. *Molecular Ecology*, 16(16), 3326-3333. https://doi.org/10.1111/j.1365-294X.2007.03394.x

Gardes, M., & Bruns, T. D. (1993). ITS primers with enhanced specificity for basidiomycetes-application to the identification of mycorrhizae and rusts. *Molecular Ecology*, 2(2), 113-118. https://doi.org/10.1111/j.1365-294X.1993.tb00005.x

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnrke, A., Rhind, N., di Palma, F., Birren, B. W., Nusbaum, C., Lindblad-Toh, K., ... Regev, A. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nature Biotechnology*, 29(7), 644-652. https://doi.org/10.1038/nbt.1883

Grigoriev, I. V., Nikitin, R., Haridas, S., Kuo, A., Ohm, R., Otillar, R., Riley, R., Salamov, A., Zhao, X., Korzeniewski, F., Smirnova, T., Nordberg, H., Dubchak, I., & Shabalov, I. (2014). MycoCosm portal: Gearing up for 1000 fungal genomes. *Nucleic Acids Research*, 42(D1), D699-704. https://doi.org/10.1093/nar/gkt1183

Grummer, J. A., Beheregaray, L. B., Bernatchez, L., Hand, B. K., Luikart, G., Narum, S. R., & Taylor, E. B. (2019). Aquatic landscape genomics and environmental effects on genetic variation. *Trends in Ecology & Evolution*, 34(7), 641-654. https://doi.org/10.1016/j.tree.2019.02.013

Harrell, F. E. Jr. (2021). *Hmisc: Harrell Miscellaneous (Version R package version 4.6-0)*. Retrieved from https://CRAN.R-project.org/package=Hmisc

Henriksen, S., & Hilmo, O. (2015). *The Norwegian red list of species*. Retrieved from Trondheim: https://www.biodiversity.no/Pages/135380/Norwegian_Red_List_for_Species

Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. *Biological Journal of the Linnean Society*, 58(3), 247-276. https://doi.org/10.1111/j.1095-8312.1996.tb01434.x

Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. *Biological Journal of the Linnean Society*, 68(1), 87-112. https://doi.org/10.1006/bijl.1999.0332

Högberg, N., Holdenrieder, O., & Stenlid, J. (1999). Population structure of the wood decay fungus *Fomitopsis pinicola*. *Heredity*, 83(3), 354. https://doi.org/10.1038/sj.hdy.6885970

Högberg, N., & Stenlid, J. (1999). Population genetics of *Fomitopsis rosea* - a wood-decay fungus of the old-growth European taiga. *Molecular Ecology*, 8(5), 703-710. https://doi.org/10.1046/j.1365-294X.1999.00561.x

Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A., & Cresko, W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. *Plos Genetics*, 6(2), e1000862. https://doi.org/10.1371/journal.pgen.1000862

Hudson, R. R., Slatkin, M., & Maddison, W. P. (1992). Estimation of levels of gene flow from DNA sequence data. *Genetics*, 132(2), 583-589. https://doi.org/10.1093/genetics/132.2.583

Mäkinen, H., Hynynen, J., Siitonen, J., & Sievänen, R. (2006). Predicting the decomposition of Scots pine, Norway spruce, and birch stems in Finland. *Ecological Applications*, 16(5), 1865-1879. https://doi.org/10.1890/1051-0761(2006)016%5B1865:PTDOSP%5D2.0.CO;2

Jombart, T., & Ahmed, I. (2011). adegenet 1.3-1: New tools for the analysis of genome-wide SNP data. *Bioinformatics*, 27(21), 3070-3071. https://doi.org/10.1093/bioinformatics/btr521

Jombart, T., Devillard, S., & Balloux, F. (2010). Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. *BMC Genetics*, 11(1), 94. https://doi.org/10.1186/1471-2156-11-94

Jones, F. C., Grabherr, M. G., Chan, Y. F., Russell, P., Mauceli, E., Johnson, J., Swofford, R., Pirun, M., Zody, M. C., White, S., Birney, E., Searle, S., Schmutz, J., Grimwood, J., Dickson, M. C., Myers, R. M., Miller, C. T., Summers, B. R., Knecht, A. K., ... Kingsley, D. M. (2012). The genomic basis of adaptive evolution in threespine sticklebacks. *Nature*, 484(7392), 55-61. https://doi.org/10.1038/nature10944

Kalliola, R. (1966). *The reduction of the area of forests in natural condition in Finland in the light of some maps based upon national forest inventories.*

Käuserud, H., & Schumacher, T. (2001). Outcrossing or inbreeding: DNA markers provide evidence for type of reproductive mode in *Phellinus nigrolimitatus* (Basidiomycota). *Mycological Research*, 105(6), 676-683. https://doi.org/10.1017/S0953756201004191

Käuserud, H., & Schumacher, T. (2002). Population structure of the endangered wood decay fungus *Phellinus nigrolimitatus* (Basidiomycota). *Canadian Journal of Botany*, 80, 597-606.

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., & Mayrose, I. (2015). CLUMPAK: A program for identifying clustering modes and packaging population structure inferences across K. *Molecular Ecology Resources*, 15, 1179-1191.

Korneliussen, T. S., Albrechtsen, A., & Nielsen, R. (2014). ANGSD: analysis of next generation sequencing data. *BMC Bioinformatics*, 15(1), 356. https://doi.org/10.1186/s12859-014-0356-4

Koskinen, M. T., Haugen, T. O., & Primmer, C. R. (2002). Contemporary fisherman life-history evolution in small salmonid populations. *Nature*, 419(6909), 826-830. https://doi.org/10.1038/nature01029

Koskinen, S., Martelin, T., Notkola, I.-L., Notkola, V., Pitkänen, K., Jalovaara, M., & Söderling, I. (2007). *Suomen väestö. (The population of Finland. In Finnish).* Helsinki: Gaudeamus.

Kuo, A., Bushnell, B., & Grigoriev, I. V. (2014). Chapter one - fungal genomics: Sequencing and annotation. In F. M. Martin (Ed.), *Advances in botanical research* (Vol. 70, pp. 1-52). Academic Press.

Lanes, É. C., Pope, N. S., Alves, R., Carvalho Filho, N. M., Giannini, T. C., Giulietti, A. M., Imperatriz-Fonseca, V. L., Monteiro, W., Oliveira, G., Silva, A. R., Siqueira, J. O., Souza-Filho, P. W., Vasconcelos, S., & Jaffé, R. (2018). Landscape genomic conservation assessment of a narrow-endemic and a widespread morning glory from Amazonian Savannas. *Frontiers in Plant Science*, 9, 532. https://doi.org/10.3389/fpls.2018.00532

Langmead, B., & Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. *Nature Methods*, 9(4), 357-359. https://doi.org/10.1038/nmeth.1923

Lilja, S., & Kuuluvainen, T. (2005). Structure of old Pinus sylvestris dominated forest stands along a geographic and human impact gradient in mid-boreal Fennoscandia. *Silva Fennica*, 39(3), 407. https://doi.org/10.14214/sf.377

Mäkinen, H., Hynynen, J., Siitonen, J., & Sievänen, R. (2006). Predicting the decomposition of Scots pine, Norway spruce, and birch stems in Finland. *Ecological Applications*, 16(5), 1865-1879. https://doi.org/10.1890/1051-0761(2006)016%5B1865:PTDOSP%5D2.0.CO;2

Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. *Cancer Research*, 27(2 Part 1), 209-220.

Maurice, S., Montes, M. S., Nielsen, B. J., Bødker, L., Martin, M. D., Jønck, C. G., Kjøller, R., & Rosendahl, S. (2019). Population genomics of an outbreak of the potato late blight pathogen, Phytophthora infestans, reveals both clonality and high genotypic diversity. *Molecular Plant Pathology*, 20(8), 1134-1146. https://doi.org/10.1111/mpp.12819
Medley, K. A., Westby, K. M., & Jenkins, D. G. (2019). Rapid local adaptation to northern winters in the invasive Asian tiger mosquito Aedes albopictus: A moving target. *Journal of Applied Ecology*, *56*(11), 2518-2527. https://doi.org/10.1111/1365-2664.13480

Meirmans, P. G., & Van Tienderen, P. H. (2004). genotype and geno-
dive: Two programs for the analysis of genetic diversity of asex-
ual organisms. *Molecular Ecology Notes*, *4*(4), 792–794. https://doi.
org/10.1111/j.1471-8286.2004.00770.x

Mela, F., Fritsche, K., de Boer, W., van Veen, J. A., de Graaff, L. H., van den Berg, M., & Leveau, J. H. J. (2011). Dual transcriptional profiling of a bacterial/fungal confrontation: Collimonas fungivorans versus Aspergillus niger. *The ISME Journal*, *5*(9), 1494–1504. https://doi.org/10.1038/ismej.2011.29

Miller, J. M., Cullingham, C. I., & Peery, R. M. (2020). The influence of a priori grouping on inference of genetic clusters: Simulation study and literature review of the DAPC method. *Heredity*, *125*(5), 269–280. https://doi.org/10.1038/s41437-020-0348-2

Murray, M. G., & Thompson, W. F. (1980). Rapid isolation of high molec-
ular weight plant DNA. *Nucleic Acids Research*, *8*(19), 4321-4325.
https://doi.org/10.1093/nar/8.19.4321

Nordén, J., Åström, J., Josefsson, T., Blumentrath, S., Ovaskainen,
O., Sverdrup-Thygeson, A., & Nordén, B. (2018). At which spa-
tial and temporal scales can fungi indicate habitat connectivity?
*Ecological Indicators*, *91*, 138–148. https://doi.org/10.1016/j.ecoli
nd.2018.03.062

Nordén, J., Penttilä, R., Siitonen, J., Tomppo, E., & Ovaskainen, O. (2013). Specialist species of wood-inhabiting fungi struggle while generalists thrive in fragmented boreal forests. *Journal of Ecology*, *101*, 701-712. https://doi.org/10.1111/1365-2745.12085

Oksanen, J., Blanchet, F., Friendly, M., Kindt, R., Legendre, P., McGlinn,
D., Minchin, P. R., O'Hara, R. B., Simpson, G. L., Solymos, P., Stevens,
M. H. H., Szoeces, E., & Wagner, H. (2019). vegan: Community Ecology
Package. 2019. R package version 2.5-6. In.

Ovaskainen, O., Schigel, D., Ali-Kovero, H., Auvinen, P., Paulin, L., Nordén, B., & Nordén, J. (2013). Combining high-throughput sequencing with fruit body surveys reveals contrasting life-history strategies in fungi. *ISME Journal*, 7(9), 1696-1709. https://doi.org/10.1038/ismej.2013.61

Penttilä, R., Lindgren, M., Miettinen, O., Rita, H., & Hanski, I. (2006). Consequences of forest fragmentation for polyporous fungi at two spatial scales. *Oikos*, *114*, 225-240. https://doi.org/10.1111/j.2006.0030-1299.14349.x

Peter, J., De Chiara, M., Friedrich, A., Yue, J.-X., Pflieger, D., Bergström,
A., Sigwalt, A., Barre, B., Freel, K., Llored, A., Cruaud, C., Labadie,
K., Aury, J.-M., Istace, B., Lebrigand, K., Barbry, P., Engelen, S.,
Lemainque, A., Wincker, P., ... Schacherer, J. (2018). Genome evo-
lution across 1,011 Saccharomyces cerevisiae isolates. *Nature,
556*(7701), 339-344. https://doi.org/10.1038/s41586-018-0030-5

Poplin, R., Ruano-Rubio, V., DePristo, M. A., Fennell, T. J., Carneiro,
M. O., Van der Auwera, G. A., Kling, D. E., Gauthier, L. D., Levy-
Moonshine, A., Roazen, D., Shakir, K., Thibault, J., Chandran, S.,
Whelan, C., Lek, M., Gabriel, S., Daly, M. J., Neale, B., MacArthur,
D. G., & Banks, E. (2018). Scaling accurate genetic variant discov-
ery to tens of thousands of samples, bioRxiv, 201178. https://doi.
org/10.1101/201178

Punttila, P., & Ihalainen, A. (2006). Luonnontilaisen kaltaiset metsät suojelu- ja ei-suojelluilla alueilla (Corrected version available from the authors). In P. Horne, T. Koskela, M. Kuusinen, A. Otsamo, & K. Syrjänen (Eds.), *METSOJät jällijlä - Etelä-Suomen metsien monimuotoisuusohjelman tutkimusraportti* (pp. 19): Ministry of Agriculture, Ministry of the Environment, Finnish Forest Research Institute and Finnish Environment Institute.

QGIS Development Team (2018). QGIS Geographic Information System
Developers Manual. Open Source Geospatial Foundation Project.

Quinlan, A. R., & Hall, I. M. (2010). BEDTools: A flexible suite of utilities for comparing genomic features. *Bioinformatics*, *26*(6), 841–842.
https://doi.org/10.1093/bioinformatics/btq033

Rochat, E., Selmoni, O., & Joost, S. (2021). Spatial areas of genotype probability: Predicting the spatial distribution of adaptive genetic variants under future climatic conditions. *Diversity and Distributions*, *27*(6), 1076–1090. https://doi.org/10.1111/ddi.13256

Rochette, N. C., Rivera-Colón, A. G., & Catchen, J. M. (2019). Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. *Molecular Ecology*, *28*(21), 4737-4754. https://doi.org/10.1111/mec.15253

Rosenberg, N. A. (2004). distract: A program for the graphical display of population structure. *Molecular Ecology Notes*, *4*(1), 137-138. https://doi.org/10.1046/j.1471-8286.2003.00566.x

Rousset, F. (1997). Genetic differentiation and estimation of gene flow from f-statistics under isolation by distance. *Genetics*, *145*(4), 1219-1228.

Rouvinen, S., Kuuluvainen, T., & Karjalainen, L. (2002). Coarse woody debris in old Pinus sylvestris dominated forests along a geographic and human impact gradient in boreal Fennoscandia. *Canadian Journal of Forest Research*, *32*(12), 2184-2200. https://doi.org/10.1139/x02-144

Savolainen, O., Lascoux, M., & Merilä, J. (2013). Ecological genomics of local adaptation. *Nature Reviews Genetics*, *14*(11), 807-820.

Schlaepfer, D. R., Braschler, B., Rusterholz, H.-P., & Baur, B. (2018). Genetic effects of anthropogenic habitat fragmentation on remnant animal and plant populations: A meta-analysis. *Ecosphere*, *9*(10), e02488. https://doi.org/10.1002/ecs2.2488

Smith, C. C., Weber, J. N., Mikheyev, A. S., Roces, F., Bollazzi, M., Kellner, K., Seal, J. N., & Mueller, U. G. (2019). Landscape genomics of an obligate mutualism: Concordant and discordant population structures between the leafcutter ant Atta texana and its two main functional symbiont types. *Molecular Ecology*, *28*(11), 2831-2845. https://doi.org/10.1111/mec.15111

Sork, V. L., Squire, K., Gugger, P. F., Steele, S. E., Levy, E. D., & Eckert, A. J. (2016). Landscape genomic analysis of candidate genes for climate adaptation in a California endemic oak. Quercus lobata. *American Journal of Botany*, *103*(1), 33-46. https://doi.org/10.3732/ajb.1500162

Stokland, J., & Kauserud, H. (2004). *Phellinus nigrolimitatus* - a wood-decomposing fungus highly influenced by forestry. *Forest Ecology and Management*, *187*, 333-343. https://doi.org/10.1016/j.foreco.2003.07.004

Storey, J. D., Bass, A. J., Dabney, A., & Robinson, D. (2019). *qvalue: Q-value estimation for false discovery rate control*. R package version 2.18.0, http://github.com/jdstorey/qvalue

Stucki, S., Orozco-terWengel, P., Forester, B. R., Duruz, S., Colli, L., Masembe, C., Negrini, R., Landguth, E., Jones, M. R., Bruford, M. W., Taberlet, P., & Joost, S. (2017). High performance computation of landscape genomic models including local indicators of spatial association. *Molecular Ecology Resources*, *17*(5), 1072-1089. https://doi.org/10.1111/1755-0998.12629

Tigano, A., & Friesen, V. L. (2016). Genomics of local adaptation with gene flow. *Molecular Ecology*, *25*(10), 2144-2164. https://doi.org/10.1111/mec.13606

Tollefsrud, M. M., Kissling, R., Gugerli, F., Johnsen, Ø., Skrøppa, T., Cheddadi, R., Van der knaap, W. O., Latałowa, M., Terhürne-berson, R., Litt, T., Geburek, T., Brochmann, C., & Sperisen, C. (2008). Genetic consequences of glacial survival and postglacial coloniza-
tion in Norway spruce: Combined analysis of mitochondrial DNA and fossil pollen. *Molecular Ecology*, *17*(18), 4134-4150. https://doi.org/10.1111/j.1365-294X.2008.03893.x

Tollefsrud, M. M., Sonstebo, J. H., Brochmann, C., Johnsen, O., Skroppa, T., & Vendramin, G. G. (2009). Combined analysis of nuclear and mi-
tochondrial markers provide new insight into the genetic structure of North European Picea abies. *Heredity*, *102*(6), 549–562. https://doi.org/10.1038/hdy.2009.16

van Boheemen, L. A., Atwater, D. Z., & Hodgins, K. A. (2019). Rapid and repeated local adaptation to climate in an invasive plant. *New Phytologist*, *222*(1), 614-627. https://doi.org/10.1111/nph.15564
Warnes, G. R., Bolker, B., Bonebakker, L., Gentleman, R., Huber, W., Liaw, A., Lumley, T., Maechler, M., Magnusson, A., Moeller, S., Schwartz, M., & Venables, B. (2020). *gplots: Various R programming tools for plotting data* (Version R package version 3.1.1). Retrieved from https://CRAN.R-project.org/package=gplots

Young, A., Boyle, T., & Brown, T. (1996). The population genetic consequences of habitat fragmentation for plants. *Trends in Ecology & Evolution, 11*(10), 413-418. https://doi.org/10.1016/0169-5347(96)10045-8

Young, A., & Clarke, G. (2000). *Genetics, demography and viability of fragmented populations* (A. G. Young & G. M. Clarke Eds.). Cambridge University Press.

Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., & Weir, B. S. (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. *Bioinformatics*, 28(24), 3326-3328. https://doi.org/10.1093/bioinformatics/bts606

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Sønstebø, J. H., Trucchi, E., Nordén, J., Skrede, I., Miettinen, O., Haridas, S., Pangilinan, J., Grigoriev, I. V., Martin, F., Kauserud, H., & Maurice, S. (2022). Population genomics of a forest fungus reveals high gene flow and climate adaptation signatures. *Molecular Ecology*, 31, 1963-1979. https://doi.org/10.1111/mec.16369