Distribution and drivers of ectomycorrhizal fungal communities across the North American Arctic

Ectomycorrhizal fungi (EMF) form symbioses with a few plant species that comprise a large fraction of the arctic vegetation. Despite their importance, the identity, abundance and distribution of EMF in the Arctic, as well as the key drivers controlling their community composition are poorly understood. In this study, we investigated the diversity and structure of EMF communities across a bioclimatic gradient spanning much of the North American Arctic. We collected roots from two principal arctic ectomycorrhizal host plants, Salix arctica and Dryas integrifolia, typically growing intermingled, at 23 locations stratified across the five bioclimatic subzones of the Arctic. DNA was extracted from ectomycorrhizal root tips and the ITS region was sequenced and phylogenetically analyzed. A total of 242 fungal Operational Taxonomic Units (OTUs) were documented, with 203 OTUs belonging to the Basidiomycota and 39 to the Ascomycota, exceeding the number of previously morphologically described EMF in the Arctic. EMF communities were dominated by a few common and species-rich families such as Thelephoraceae, Inocybaceae, Sebacinaceae, Cortinariaceae, and Pyronemataceae. Both host plants showed similar species richness, with 176 OTUs on Salix arctica and 154 OTUs on Dryas integrifolia. Host plant identity did not affect EMF community composition. The ten most abundant OTUs had a wide geographic distribution throughout the Arctic, and were also found in boreal, temperate and Mediterranean regions, where they were associated with a variety of hosts. Species richness did not decline with increasing latitude. However, EMF community structure changed gradually across the bioclimatic gradient with the greatest similarity between neighboring bioclimatic subzones and locations. EMF community structure was correlated with environmental factors at a regional scale, corresponding to a complex of glaciation history, geology, soil properties, plant productivity and climate. This is the first large-scale study of EMF communities across all five bioclimatic subzones of the North American Arctic, accompanied by an extensive set of environmental factors analyzed to date. While our study provides baseline data to assess shifts of plant and fungi distribution in response to climate change, it also suggests that with ongoing climate warming, EMF community composition may be affected by northward shifts of some taxa.


INTRODUCTION
Ectomycorrhizal fungi (EMF) are critical to growth and survival of their host plants because they provide water and limiting nutrients in exchange for photosynthetic carbon (Smith and Read 2002).They are thought to be particularly important in harsh environments such as the Arctic tundra, a permafrost-underlain cold-dominated biome with low nutrient availability north of the treeline.EMF have been studied up to 798 N where they were found to co-occur with their host plants (Bledsoe et al. 1990, Kohn andStasovski 1990).Although in the Arctic EMF associate only with 6% of the vascular plant species, including shrubs, a few sedges and forbs (e.g., Dryas, Salix, Betula and Polygonum), these plants are important.For example, they are dominant species in plant community types that cover up to 69% of the ice-free Arctic (Walker et al. 2005).
The Arctic is divided into five bioclimatic subzones (A-E), with A being the coldest and E being the warmest (Walker et al. 2005).With increasing latitude, organisms face harsher conditions due to decreasing air and soil temperatures (Billings 1992).Vascular plant species richness, including that of EMF host plants, decreases.Simultaneously, plant communities change from a zone of low shrubs in the Low Arctic (subzone E) to a zone of cushion forbs without shrubs in the High Arctic (subzone A) (Walker 2000).This change in plant community composition across the bioclimatic gradient is driven primarily by climate and soil pH (Walker et al. 2011).
Recent climatic changes in the Arctic have led to a pan-arctic shrub expansion (Tape et al. 2006) and an increase in tundra productivity and greening of the Arctic (Bhatt et al. 2010).With continuing warming, further large changes in arctic plant communities are expected (Callaghan et al. 2004).Paleobotanical studies show that plant and fungal communities underwent major changes during past glacial and interglacial cycles, with an increase in shrubs and trees and their root associated EMF since the last glaciation (Lydolph et al. 2005, de Vernal andHillaire-Marcel 2008).Experimental warming of low arctic tundra has shown not only a shift of plant communities toward increased shrub dominance (Clemmensen and Michelsen 2006, Walker et al. 2006, Elmendorf et al. 2012), but also changes in the community structure of their associated EMF (Deslippe et al. 2011).EMF are expected to play an important role in these ongoing climatedriven changes in plant communities, particularly shrub expansion, since Arctic shrub species are all ectomycorrhizal.
To understand the factors underlying shrub expansion in the Arctic, we need to include studies of their EMF.In particular, we need to know their identity, their distribution and what factors shape these communities in the Arctic.However, despite the ubiquity and significance of EMF for Arctic tundra and the fact that they have been reported for more than a century (Hesselman 1900, Bledsoe et al. 1990, Kohn and Stasovski 1990, Gardes and Dahlberg 1996, Newsham et al. 2009), our knowledge of the identity, distribution and ecology of arctic EMF is only fragmentary (reviewed in Timling and Taylor 2012).
Thus far, arctic EMF studies have relied primarily on sporocarp (mushroom) collections of macrofungi.Additionally, studies have been limited in geographic and temporal scope because the Arctic is not easily accessible and occurrence of sporocarps is ephemeral and irregular.Lists of macrofungi from Greenland, Svalbard, the Russian Arctic and Iceland report about 2600 morphologically described macrofungi, with at least 150 ectomycorrhizal species (Elvebakk and Prestrud 1996, Karatygin et al. 1999, Hallgrimsson and Eyjolfsdottir 2004, Borgen et al. 2006).Such lists have not yet been compiled for North America.Patterns from morphological sporocarp-descriptions of macrofungi indicate that arctic fungi are widely distributed in arctic and alpine habitats on all continents.Some widely distributed EMF genera that have a preponderance in arctic and alpine conditions include Inocybe, Cortinarius, Hebeloma, Russula, Thelephora, Tomentella, Cenococcum, and Laccaria (Gardes and Dahlberg 1996, Mu ¨hlmann and Peintner 2008, Ryberg et al. 2009, Deslippe et al. 2011, Fujiyoshi et al. 2011).However, reliance on morphologically recognized species may underestimate fungal diversity (Taylor et al. 2006).Hence, the reported species richness for arctic EMF provides only a conservative estimation.
Outside the Arctic, many fungi have distribution patterns corresponding with geography (e.g., Taylor et al. 2006, Geml et al. 2008b, O'Donnell et al. 2011).Distinctive phylogeographic patterns have been observed for fungi in boreal, temperate and tropical regions (Taylor et al. 2006, Geml et al. 2008a).In the Arctic however, studies of lichen-mycobionts and EMF have not revealed discrete distributions at the continental or smaller scales, but instead highlight the importance of wide-ranging dispersal (Buschbom 2007, Geml et al. 2010a, Geml et al. 2012b).Thus we might expect that species of EMF in the Arctic are widely distributed.A wide distribution of EMF would be facilitated by little or no host preference.Indeed low host preference has been shown for EMF communities associated with several alpine and arctic dwarf shrubs (Kernaghan and Harper 2001, Ryberg et al. 2009, Ryberg et al. 2010, Fujimura and Egger 2012).However, most of these studies were conducted on local (;5 km) scales, and hence require confirmation at a much broader geographic scale.In contrast to the relatively low diversity of EMF host plants in the Arctic, studies exploring EMF diversity in the Arctic using molecular methods (Clemmensen and Michelsen 2006, Wallenstein et al. 2007, Fujimura et al. 2008, Bjorbaekmo et al. 2010, Deslippe et al. 2011, Fujiyoshi et al. 2011, Geml et al. 2012b) have revealed species-rich EMF communities; this parallels findings from temperate and boreal biomes (Tedersoo et al. 2006, Buee et al. 2009, Taylor et al. 2010).The largest-scale study to date, exploring EMF diversity of D. octopetala across two of the five bioclimatic subzones in the European Arctic, observed no decline in species richness with increasing latitude (Bjorbaekmo et al. 2010).Whether this pattern applies across the five bioclimatic subzones and to other dwarf shrubs and regions of the Arctic has yet to be determined.
The objectives of our study were (1) to characterize patterns of EMF community composition along a bioclimatic gradient across all five subzones of the North American Arctic and (2) to identify key drivers controlling their community composition.We sampled two widespread and important shrubs of the arctic tundra, Salix arctica Pall and Dryas integrifolia Vahl (Hulte ´n 1968) in the North American Arctic and carried out detailed molecular phylogenetic analyses of EMF colonizing their root tips.
To identify patterns of EMF community composition along the bioclimatic gradient in the Arctic we hypothesized the following: (1) EMF richness decreases with latitude, parallel to species richness of macro-organisms.(2) EMF communities of S. arctica and D. integrifolia change gradually across bioclimatic subzones, as observed in plant communities.(3) Dominant EMF associates of S. arctica and D. integrifolia have wide distributions across Arctic and Alpine habitats.However, these dominants will be restricted to cold-dominated regions due to cold adaptation.(4) Climate and abiotic soil factors are the most important drivers of EMF community composition, while host plant identity of S. arctica and D. integrifolia exerts little influence.In our study we have sampled for the first time all five bioclimatic subzones of the Arctic and analyzed key drivers of EMF community structure based on a wide range of environmental factors along the gradient.

MATERIAL AND METHODS
This paper reports the combined work of two research initiatives conducting descriptive surveys in the North American Arctic.Because the collaboration was established after sampling and molecular analysis had been completed, we provide information regarding these steps separately, where they differ.The first dataset comprises sampling sites throughout the Canadian Arctic Archipelago (CAA) (Eriksen et al. 2006), and the second dataset consists of sampling sites along the North American Arctic Transect (NAAT) (Walker et al. 2008).
In both datasets we sampled sites in five bioclimatic subzones of the Arctic (A-E) as portrayed on the Circumpolar Arctic Vegetation Map (Walker et al. 2005) (Fig. 1).A 'bioclimatic subzone' is defined based on a combination of summer air temperature and dominant plant growth forms.Subzone A is the coldest, with a Mean July Temperature (MJT) of 0-38C.The dominant plant growth forms are cushion forbs, mosses and lichens.This coldest subzone has a very limited extent.Subzone B has a MJT of 3-58C and is characterized by prostrate dwarf shrubs.Subzone C has a MJT of 5-78C, and is dominated by hemi-prostrate dwarf shrubs, sedges and mosses.Subzone C represents the subzone with the largest extent in the circumpolar Arctic tundra.Subzone D has a MJT of 7-98C, with erect dwarf shrubs, tussock sedges and mosses.Subzone E is the most southern subzone and has a MJT of 9-128C, with low shrubs, tussock sedges and mosses.Vascular plant diversity and plant cover increases from subzone A to E (Walker et al. 2005).
To analyze the data, we structured our observations into two datasets.The first consisted of all data collected along both the CAA and the NAAT.We refer to this as the ARCTIC dataset.This dataset represents a large geographical scale, but is described with fewer environmental factors than the second dataset, which consists of data collected only along the NAAT.While this dataset contains fewer sites in each bioclimatic subzone, it is accompanied by a more extensive collection of environmental data for each study site.
Due to the different strengths of these two datasets, we used the ARCTIC dataset to assess the overall EMF diversity associated with S. arctica and D. integrifolia, the species richness along a latitudinal gradient, the geographical distribution patterns of the observed fungal taxa, and the effect of host plant and environment on the ECM fungal communities.We used the NAAT dataset to investigate changes in EMF community structures between sampling locations and bioclimatic subzones.Finally, we used the NAAT dataset to identify key environmental factors affecting fungal communities along this bioclimatic gradient.

Study area
NAAT.-We sampled ectomycorrhizal root tips associated with S. arctica and D. integrifolia at eight sites along the 1800 km North American Arctic Transect (NAAT), covering bioclimatic subzones A-E (Table1).All sites represent mesic zonal sites that were extensively studied by the 'Biocomplexity Project' (Walker et al. 2008).Detailed descriptions of the study sites are provided in Kade et al. (2005), Vonlanthen et al. (2008), andRaynolds et al. (2008).Samples along the NAAT were collected at the end of July in 2005, 2006 and 2009.In 2006, we obtained samples of D. integrifolia and S. arctica from Thule (Greenland).This location was a mesic site from subzone C, and can be described as prostrate dwarf shrub herb tundra.

Sampling and processing
CAA.-The material for this study was collected July 1st to September 1st 1999 as part of the Tundra North West 1999 expedition in the Canadian Arctic, spanning 14 sites (Eriksen et al. 2006).These sites were selected to represent longitudinal and latitudinal gradients, to encompass vegetation spanning the Low to High Arctic and four bioclimatic subzones (B-E).Sites included mesic and dry conditions (Table 1).Further information can be found in Eriksen et al. (2006).
CAA.-At five sites, two intensively studied 20 3 20 m plots were established, with one in mesic conditions and one in dry conditions.These plots were not more than 500 m apart.Each plot was surveyed for vegetation and soil characteristics within the Biodiversity Program of Tundra North West 1999 (Bo ¨lter 2006, Bo ¨lter et al. 2006, Eriksen et al. 2006), and eight plants of both S. arctica and D. integrifolia, including their full root systems were randomly collected from each plot (Table 1).In order to obtain additional mycorrhizal samples from less intensively studied sites and to cover a larger area and hence a larger potential proportion of local species richness, 2-10 plants of S. arctica were randomly collected along 500 m transects at each of the fourteen sites and 3-10 plants of D. integrifolia were collected at three sites along the same transect.At the five intensively studied sites, the transects were an outward extension beyond the plot area.In total, 164 S. arctica and 89 D. integrifolia root systems were obtained.Roots were processed in a similar fashion as samples from the NAAT.
Single or multiple representative root tips of all distinguished morphotypes from each plant were sorted and placed individually in 1.5 ml tubes.From each plant, at least 5 mycorrhizal tips were collected, regardless of whether one or several morphotypes were detected.All discernible mycorrhizal morphotypes were sampled from each plant.Immediately after detaching, each mycorrhizal root tip was submerged in 300 lL of 2X CTAB lysis and kept frozen at À208C until DNA extraction was performed as described by Gardes and Bruns (1996).
NAAT.-At each site, we randomly chose three to eight plants from each shrub species in the immediate vicinity of the 10 3 10 m grid established for the 'Biocomplexity of Frost-Boil Ecosystems Project', representing mesic zonal conditions (Kade et al. 2005, Raynolds et al. 2008, Vonlanthen et al. 2008, Walker et al. 2008).Entire plants were excavated along with soil surrounding the root system, for a total of 34 S. arctica and 39 D. integrifolia (Table 1).The plants were stored in a cooler and transported to the laboratory within ten days, where the root systems were immediately washed and rinsed with de-ionized water.All root tips were harvested and morphotyped following Agerer (1987Agerer ( -2002)).Each morphotype was kept separate for each plant throughout analysis (i.e., no attempt was made to define shared morphotypes across samples because this is an error-prone practice (Sakakibara et al. 2002).Two root tips from each morphotype were placed individually in cryovials and frozen at À808C, then lyophilized prior to DNA extraction.

Molecular analysis
CAA: DNA extraction and PCR.-Cell disruption was performed by three freeze/thaw cycles in liquid nitrogen/658C.EMF root tips were crushed with a motorized mini-pestle (Kontes, Vineland, NJ, USA) before incubation at 658C for 1 h.After chloroform extraction and DNA precipitation of 1856 samples, the entire ITS region was amplified using fungal specific primer ITS1-F (Gardes and Bruns 1993) and generic primer ITS4 (White et al. 1990) with a PTC-200 thermal cycler (MJ Research, Waltham, MA, USA) following the protocol by Gardes and Bruns (1993).
NAAT: DNA extraction and PCR.-The lyophilized roots were ground with 3.2 mm stainless steel beads (BioSpec.Products) using a DNA Mixer Mill (Retsch, Haan, Germany).Ground root tips were immersed in buffer and RNase for 12 h before DNA extraction to increase yield.DNA was extracted from 776 single root tips, using the E.Z.N.A. SP Fungal DNA Kit (Omega Bio-tek, Doraville, GA).Extracted DNA was amplified with primers ITS1-F and ITS4 as above using a PTC-225 Thermal Cycler (MJ Research, Waltham, MA, USA).The polymerase chain reaction (PCR) program was 2 min at 968C, then 34 cycles of 30 s at 948C, 40 s at 548C, 1 min at 728C followed by 10 min at 728C.

RFLP screening
PCR products (700-800 bp) were digested with the restriction enzyme HaeIII (BioLabs, New England, USA) following manufacturers instructions and visualized on the 2% Nusieve GTG plus 1% agarose gels.Restriction fragment length polymorphism (RFLP) patterns were identified using AlphaImager 3400 (Alpha Innotech, San Leandro, CA, USA).Only bright bands that added up to the correct total fragment size for the ITS region were included in the designation of particular RFLP patterns, while weak bands were excluded.RFLP profiles were grouped with the program GERM 1.01 (Dickie et al. 2003).To minimize the underestimation of EMF-diversity, we tabulated RFLP types for each plant species separately and matched RFLP types only within a site.We applied the following settings: maximum forward and backward error 10 bp, maximum sum error 100 bp, lower measurement limit 100 bp, and the joint matching and ranking method (Dickie et al. 2003).For the CAA mycorrhizae, 71% of the sequenced samples were digested with four restriction enzymes (Cfo1, Hinf1, HaeIII, and MboI) independently, 24% with two enzymes (CfoI, HinfI), and 4% with three enzymes (CfoI, HinfI, HaeIII).Identical RFLPpatterns were then used to molecularly characterize each mycorrhizae.
We recorded the number and abundance of each RFLP type for each plant at each site.Multiple samples from different sites for each RFLP type were chosen for sequencing.More sequences were obtained for frequent RFLPtypes, up to a maximum of 22 sequences from 10 sites for one RFLP-type.Altogether, this allowed us to later combine RFLP types from different sites and plant species having similar or identical sequences into a single operational taxonomic unit (OTU).

Cloning and sequencing
RFLP types indicating colonization of a root tip by a single fungal taxon were purified using QIAquick PCR Purification Kit (Qiagen Sciences, Valencia, CA, USA).Purified DNAs were cycle sequenced with Applied Biosystems (ABI) Big Dye v.3.1 kit using ITS1-F and ITS4.Cycle sequence products were purified using Sephadex-G50 (GE Healthcare-Bio Sciences-AB, Uppsala, Sweden) before they were sequenced on a capillary DNA sequencer 3130x/Genetic Analyzer (Applied Biosystems, USA) at the University of Alaska (Fairbanks, USA) or Genome Express (Meylan, France).
NAAT samples that contained multiple taxa were cloned as follows.PCR products were purified and concentrated with Zymo DNA Clean and Concentrator-5 columns (Zymo Research, Irvine, CA, USA).DNA concentration was determined with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Rockland, DE, USA).For cloning we used the TOPO TA Cloning Kit for Sequencing (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions.From each sample/plate, five clones were randomly chosen and PCR amplified with M13F and M13R primers.PCR products were digested with HaeIII.From these five clones, one representative for each unique clone RFLP type was sequenced as above using the primers M13F and M13R.All sequences were deposited in GenBank under the accession numbers JX630331-JX630750 and JX630751-JX630968.

Bioinformatics and statistical analysis
Sequence quality filtering.-Sequenceswere cleaned according to Taylor and Houston (2011).In brief, vector sequences were removed and the two reads for each sequence were assembled using Aligner v.1.3.4 (Codoncode, Dedham, MA, USA).Consensus sequences were exported in FASTA format.Base calls with a phred score ,20 were masked with the letter 'N' using a Perl script.Ambiguous ends of sequences were trimmed using the program TRIMSEQ (EMBOSS) at http://cbi.labri.fr/utils/Pise/trimseq. html, with a window size set to 20 and a threshold ambiguity set to 5%.Trimmed sequences were submitted to the program Bioinformatics ToolBox (DNA 2.0) at http://www.dna20.com/index.php?pageID¼151 in order to identify and eliminate sequences with .2%Ns.

OTU clustering and alignment
In order to assess phylogenetic diversity based on the ITS1-5.8S-ITS2region, we clustered the sequences into OTUs using CAP3 (Huang and Madan 1999).OTUs were defined as sequences sharing 97% sequence identity in the overlapping regions.We set 'overlap percent identity' to 97, 'maximum overhang percent length' to 60, and 'clipping range' to 6.All other settings were left at defaults.For each OTU, the number and the sampling location of clones and matching RFLP patterns were recorded.For further community analysis, we recorded the presence/absence of an OTU in each plant replicate at a sampling site.For example, if 8 Salix plants were sampled from a particular site, the range of possible abundances for an OTU was 0 to 8. We included all sequences that fell within an OTU for phylogenetic analyses.
In order to identify the OTUs, we compared the ITS sequences to a curated database derived from GenBank sequences hosted on the Fungal Metagenomics website (http://biotech.inbre.alaska.edu/fungal-metagenomics/),using the FASTA matching algorithm (Pearson 2000).From these search results, named sequences with the highest similarity were included in sequence alignments and phylogenetic analyses.Sequences were considered intergeneric chimeras when they had .90%sequence similarity to different genera in ITS1 and ITS2 and intrageneric chimeras when ITS1 and ITS2 had .97%sequence similarity to different species.
Because the ITS regions of different genera cannot be aligned, we grouped the sequences based on their BLAST search result into genera/ families.OTU sequences of the various genera were aligned in Clustalw (Li 2003).Initially we built a maximum parsimony tree in PAUP* 4b10 (Swofford 2003) to identify major clades within genera that were too divergent to be aligned with confidence.Sequences from the resulting major clades were aligned separately in Clustalw, as above.The alignments were further adjusted by eye in the alignment editor Se-Al (Rambaut 1996).PAUP* 4b10 (Swofford 2003) and Modeltest 3.7 (Posada and Crandall 1998) were used to determine the best fit evolutionary model, with the Hierarchical Likelihood Ratio Test (hLRT).Trees were constructed using the maximum likelihood method in Garli 0.951 (Zwickl 2006).Statistical branch support was estimated using 100 maximum likelihood bootstrap replicates.A consensus tree was computed based on 50% majority rule.

Nomenclature
OTUs were named on the basis of the phylogenetic results.OTUs were considered to be identified to species level when the closest match in the 90% bootstrap supported clade was a single taxon with .97%sequence similarity (i.e., all clade members represent the same species, or include one fully identified species and other unidentified taxa).We applied 'affinity' (aff.) when the closest matches were fully identified (genus þ species) and showed 93-97% sequence similarity and there was no species incongruence among the members of the clade in regard to the species name.We used 'species' (sp.) when the closest matches were only described at the genus level, as well as when the closest matches with .95%similarity were fully identified, but had different species names (incongruence in the terminal clade).
A sequence similarity below 93% was used to assign taxa at the family level, and below 83% at the order level (Deslippe et al. 2011).

Ordination analysis
We used nonmetric multidimensional scaling (NMDS) to investigate relationships between the distributions of ectomycorrhizal fungi and envi-ronmental factors.Analyses were executed in PC-ORD 5 (McCune and Grace 2002).We used the abundance-based version of the Sorensen index (Bray-Curtis) as a distance measure.Prior to ordination, we conducted outlier tests for OTUs.If the deviation exceeded 2 SD, an OTU was considered an outlier and was excluded from the analysis.We determined the dimensionality of the ordination and chose the lowest dimensionality that captured most of the variation.Both the ARCTIC and NAAT datasets were best described by 3-dimensional solutions that had instabilities below 0.00001.To avoid solutions involving local minima in stress values, each analysis was run 50 times using a random seed, 50 runs of real data, and 500 iterations with Monte Carlo randomizations to test for significance.
The environmental factors included in the analysis consisted of measurements taken at the ground, as well as remotely-sensed data from satellites.Categorical factors included 'host plant' and 'location'.Numerical environmental factors are listed in Table 3.Before NMDS of the ARCTIC dataset, we square-root transformed the main matrix to down-weight the importance of abundant OTUs (McCune et al. 2002).In order to investigate whether ectomycorrhizal fungal communities differed according to location (23 locations), bioclimatic subzone (A-E), host plant (S.arctica, D. integrifolia), or soil moisture (mesic, dry), we used non-parametric multiple response permutation procedures (MRPP) in PC-ORD 5. MRPP determines whether defined groups vary statistically from those derived by random assembly.It measures within-group homogeneity, A (analogous to 'effect size'), which can reach a maximum of 1, when all communities within a group are identical and have no overlap with other groups.In contrast, if within-group variation occurs by random chance, then A ¼ 0. If A .0, there is among-group variation that is not completely explained by within-group variation, and we can calculate the probability that group differences are due to chance (McCune et al. 2002).We report the chance-corrected effect size.Because we used MRPP in conjunction with NMDS ordination, we chose the same distance measure Sorensen (Bray Curtis), as recommended by McCune et al. (2002).
To illustrate the changes of EMF community composition across the bioclimatic subzones, we applied a two-way cluster analysis in PC-ORD 5.The two-way cluster dendrograms present OTUs of selected genera that occurred more than once in the ARCTIC across the five bioclimatic subzones.

Diversity analyses
Due to the different number of samples taken at various sampling sites, we randomly subsampled three plants to estimate measures of species richness.Species richness (Mao Tau, Chao2) of the fungal communities were calculated in EstimateS 8.0 (Colwell 2006).Rarefaction curves were also computed for each bioclimatic subzone using EstimateS 8.0.The number of OTUs was determined by randomly subsampling the observed OTUs 50 times.
To examine the relationship between the observed and estimated OTU richness (Mao Tau, Chao2) along the latitudinal gradient, we performed a linear regression analysis in R (R Development Core Team 2008).

EMF diversity in the Arctic
To identify EMF diversity in the Arctic, we sampled 326 plants (198 S. arctica and 128 D. integrifolia) across the North American Arctic.Following morphotyping of root tips, extraction of DNA and RFLP analysis we obtained 644 fungal, non-chimeric sequences.Clustering these sequences across the entire ITS region resulted in 242 OTUs with 203 OTUs belonging to the Basidiomycota and 39 to the Ascomycota.When considering the abundance of OTUs (based on the presence/absence of each OTU across all plants), we found that the 203 Basidiomycete OTUs occurred 751 times, whereas the 39 Ascomycete OTUs occurred 91 times (Supplement).The rank abundance curve shows a lognormal distribution, with a few abundant OTUs and many rare species.In fact, 111 OTUs (46%) were singletons (Appendix: Fig. A1).We identified all OTUs in 37 separate maximum likelihood phylograms (Fig. 1-37) at http://dx.doi.org/10.5061/dryad.ff1g6.
Overall biodiversity was dominated by only a few fungal families and orders, which included Thelephoraceae, Inocybaceae, Sebacinaceae, Cortinariaceae and Pyronemataceae (Fig. 2).Twenty nine OTUs (12%) were identified to the species level, 194 OTUs (80%) to the genus level and 19 OTUs (8%) could be identified to family or order.Fifty eight OTUs (24% of all OTUs) showed ,95% similarity to any publicly available sequences on GenBank (as of December 2010).Of these, 44 OTUs were Basidiomycota and 14 were Ascomycota.The genus Inocybe had the most (21 OTUs) sequences without matching reference sequences in GenBank.
The ten most abundant OTUs (number of OTUs in parentheses) belonged to the genera Cortinarius (5), Hebeloma (1), Tomentella (1), Sebacina (1), Inocybe (1) and Laccaria (1).Their phylogenetic relationships to the closest matches in and outside the Arctic are shown at http://dx.doi.org/10.5061/dryad.ff1g6(their Figs.1,7,9,12,13,14,21,23).They were found in 7-16 of the 23 sampling sites.The geographic distributions of four of the ten most abundant OTUs are shown in Fig. 3. BLAST-searches in GenBank revealed that all of the most abundant taxa had close matches outside of the Arctic, and they also occurred at other Arctic locations.Additionally, these OTUs were found in a variety of habitats from the Arctic to the Mediterranean and associated with a range of host plants, including angio-and gymnosperm tree species, orchids, a sedge and a perennial forb (Table 2).Overall, 73% of all OTUs (176) had matches in GenBank with .97%similarity in the ITS region that originated in other regions of the Arctic and beyond the Arctic.
The six OTUs found on S. arctica in subzone A included Tomentella sp.27, T. stuposa 1, T. aff.terrestris (OTU 22,25,188), Sebacina sp.23 (OTU 21), Cortinarius sp. 15 (OTU 39) and Laccaria sp.(OTU 20).These OTUs occurred from one subzone to all five subzones and also on D. integrifolia.Closest matches (!97% similarity) in GenBank showed that most of them occurred outside the Arctic on a variety of host plants.In fact, OTUs 39 and 20 were among the ten most abundant OTUs throughout our study area, and OTU 39 in particular had a wide geographic distribution outside the Arctic and occurred on a wide array of host plants (Table 2).These results indicate that neither ectomycorrhizal fungi in High or Low Arctic are limited to the Arctic, to arctic/alpine (tundra) habitats, nor restricted to a particular host (Salix, Dryas).

Patterns of EMF communities
Species richness along the latitudinal gradient.-Linearregression analysis of observed and estimated species richness (Mau Tao, Chao2) showed a non-significant decrease in EMF species richness (OTUs) associated with S. arctica and D. integrifolia with increasing latitude (Fig. 4).
Community structure across the bioclimatic subzones.-NMDSordination of the NAAT dataset revealed changes in EMF community structure across the bioclimatic subzones (Fig. 5A).Additionally, two way cluster analysis supports the observation that neighboring subzones are most similar.Only OTU 39 (Cortinarius sp.15) occurred across all bioclimatic subzones while many OTUs (belonging to the genera Tomentella, Thelephora, Cortinarius, Inocybe, Entoloma, Russula, Lactarius, Clavulina, Sebacina, Cenococcum) were just observed in subzone C.However, the indicator species test did not show any species that were uniquely associated with particular bioclimatic subzones.This result might be due to low sampling intensity in our study and therefore lack of statistical power.
Despite the lack of indicator species, MRPP analysis showed a significant effect of 'bioclimatic subzone' on EMF community structure along the NAAT (A ¼ 0.037, P , 0.001) and across the ARCTIC dataset (A ¼ 0.007, P , 0.001).To refine these differences in community structures, we considered the distribution of OTUs from the most abundant genera, across the subzones B-E (excluding subzone A, because it was represented by only three plants of S. arctica).Using MRPP we found that OTUs belonging to Tomentella were significantly different among all four subzones (A ¼ 0.032, P , 0.001; Appendix: Fig. A2A).Sebacina had no differences in the Low Arctic subzones (D-E), while differences occurred between the High Arctic subzones (B-C) and between the High and Low Arctic subzones (A ¼ 0.026, P , 0.001) (Appendix: Fig. A2B).Inocybe showed mixed results, with no difference between the High Arctic subzones (B-C), and also no differences between the Low Arctic subzone E with the High Arctic subzones B and C. Nevertheless, Inocybe differed among subzones in the Low Arctic (A ¼ 0.015, P , 0.001; Appendix: Fig. A2C).In contrast, OTUs of Cortinarius did not differ among the subzones (A ¼ 0.001, P , 0.556; Appendix: Fig. A2D).In summary, across the bioclimatic subzones, we find opposing distribution patterns between Tomentella and Cortinarius (different community structures among all subzones for Tomentella and none for Cortinarius) and between Sebacina and Inocybe (different community structures between the High and Low Arctic subzones for Sebacina and none for Inocybe).
Environmental drivers.-Three-dimensionalNMDS ordinations were used to describe variation in EMF community structure across the Arctic.The ARCTIC dataset consisted of 239 OTUs spanning 38 sites (23 locations 3 2 host plants), after three outliers were removed.Five sites were only represented by one host plant species.The three axes accounted for 65.7% of the variation in community composition.Axis 1 was correlated with temperature-related factors such as subzone and Summer Warmth Index (SWI), which is the sum of the mean monthly temperatures above 08C (Walker et al. 2011).Axes 2 and 3 correlated with SWI, geographic location ('longitude'), and geology-related factors ('landscape age', 'pH') (Fig. 5B).
The NAAT dataset consisted of 128 OTUs from 12 sites (6 locations 3 2 host plants), after removing seven outliers.The environmental matrix included 35 numerical (Table 3) and two categorical factors.A three dimensional NMDS ordination was used to describe variation in fungal community structure along the NAAT (Fig. 5A).The three axes accounted for 70.1% of the variation in the dataset.Axis 1 correlated to soil moisture, while axis 2 corresponded with vegetation factors, including an increase in vascular plant species richness and plant pro-ductivity.Other contributors to axis 2 included soil carbon content, and summer and winter precipitation to the south along the NAAT.Axis 3 was best interpreted as a complex of factors relating to geographic location, soil chemistry and air and soil temperature (Table 3).MRPP analysis of dry and mesic sites across five locations (MI, BI2, EL, DI, BFS) showed that EMF communities of S. arctica and D. integrifolia were weakly but significantly (A ¼ 0.005, P ¼ 0.001) affected by these two categories of soil moisture.The NMDS ordinations of both datasets (ARCTIC, NAAT) showed that EMF communities along the bioclimatic gradient in the Arctic shift along the bioclimatic gradient, mainly in concert with a complex set of temperature and soil factors associated with each geographic location.

EMF diversity in the Arctic
We found diverse EMF communities associated with S. arctica and D. integrifolia in the North American Arctic.The observed 242 OTUs on two hosts exceed the number of previously reported sporocarps and morphologically described mycorrhizas (approximately 150) in the Arctic across all regions and host species.With 176 OTUs on S. arctica and 154 OTUs on D. integrifolia, we found a species (OTU) richness similar to that previously observed on D. octopetala (137 OTUs) along a latitudinal gradient from Alpine Norway to Svalbard (subzones B and C) (Bjorbaekmo et al. 2010).As in most microbial studies, our rarefaction curves did not approach an asymptote and 46% of the OTUs were singletons.This result indicates that we did not exhaust EMF species richness and that the actual diversity is certainly greater than observed in our study.
The major genera found in our study (Thelephora, Tomentella, Sebacina, Inocybe, Cortinarius, Russula, Hebeloma, Laccaria, Clavulina) are characteristic of arctic and alpine environments (Mu ¨hlmann and Peintner 2008, Ryberg et al. 2009, Bjorbaekmo et al. 2010, Deslippe et al. 2011, Geml et al. 2012b).While the Thelephoraceae have a world-wide distribution (Koljalg et al. 2000), they seem especially species-rich in the Arctic, as they comprised nearly a third of all observed species in our study.Moreover, the proportional increase of melanized fungi (e.g., Thelephoraceae) from warm to cold subzones in our study corresponds with the increase of thelephoralean fungi from central Norway to Svalbard (Bjorbaekmo et al. 2010), demonstrating the fitness of this fungal group to the arctic environment.
In contrast to other studies in the Arctic (Hrynkiewicz et al. 2009, Bjorbaekmo et al. 2010, Fujiyoshi et al. 2011), Cenococcum geophilum and species in the genera Russula and Lactarius were less abundant with latitude along the bioclimatic gradient.We found Cenococcum geophilum only in subzone C; it is often associated with high carbon content in the soil, which in fact was highest in subzone C in our study.The genera Russula and Lactarius are abundant and species-rich EMF in the Low Arctic and the Boreal Forest (Geml et al. 2009, Geml et al. 2010b, Taylor et al. 2010, Deslippe et al. 2011), while they were less abundant and species-rich in our higher latitude study.This suggests, that these genera may prefer EMF host plants of the Low Arctic and Boreal and/or may be less adapted to the climatic conditions in the colder subzones of the Arctic.
Notably, members of the Russulaceae, Thelephoraceae, and Cenococcum are ubiquitous and dominant components of ECM communities in temperate forests throughout the world (e.g., Gardes andBruns 1996, Jonsson et al. 1999, Table 3. Correlation coefficients (r 2 ) for variables in the NMDS ordination along the NAAT and throughout the ARCTIC.Normalized Difference Vegetation Index (NDVI) is an index of vegetation greenness and is commonly used as indicator for biomass (Walker et al. 2011).

Group
à Summer Warmth Index (SWI) is the sum of mean monthly temperature .08C(Walker et al. 2011).§ Thawing Degree Days (TDD) is the sum of mean daily temperatures .08Cover a year (Walker et al. 2011).} Freezing Degree Days (FDD) is the sum of mean daily temperature ,08C over a year (Walker et al. 2011).Tedersoo et al. 2006).Species of Inocybe and Cortinarius do occur in temperate regions, but do not appear to be ubiquitous dominants in warmtemperate regions, while they become increasingly common in the boreal forests of Europe and North America (Jonsson et al. 1999).These patterns suggest that, at the genus level, Tomentella-Thelephora are climate generalists, while Inocybe and Cortinarius may be more specialized to cold climates and the Russulaceae to warmer climates, though extending to boreal regions.

Large scale EMF community patterns in the Arctic
Wide distribution of arctic EMF.-The majority (73%) of OTUs in our study matched GenBank sequences recovered from other regions both within and beyond the Arctic.Thus, our study supports patterns from sporocarp records that have claimed a wide distribution of arctic fungi across arctic and alpine habitats spanning all continents (Ronikier and Ronikier 2010).It also agrees with molecular studies showing a panarctic distribution of Lichenomphalia (Geml et al. 2012a) and another study, which showed that 73.3% of ectomycorrhizal phylotypes observed on the isolated Arctic archipelago of Svalbard also occurred outside of Svalbard (Geml et al. 2012b).
Analysis of the ten most abundant OTUs showed that they were not restricted to arctic and alpine regions.Instead they also occurred in boreal, temperate and Mediterranean regions, in a wide variety of habitats and with a wide range of host-plants (Table 2).This suggests either that these fungal species have very wide ecological amplitudes and niches or that there has been recent population differentiation which is not reflected among our ITS groupings at 97% similarity.
The distribution patterns we observed suggest the contribution of both terrestrial and transoceanic dispersal over long distances across multiple scales to the assembly of arctic communities.While individual long-distance dispersal events (especially transoceanic) are considered rare, it has occasionally been hypothesized for hypoand epigeous fungi at lower latitudes (e.g., Halling et al. 2008, Hosaka et al. 2008).Sufficient dispersal over large distances to homogenize populations has been demonstrated to occur in some lichens (Buschbom 2007, Geml et al. 2010a) and plants occurring in the Arctic (Tremblay andSchoen 1999, Alsos et al. 2007).Vectors for such dispersal of plants over wide areas of the Arctic include wind, snow, birds, driftwood, sea ice and mammals (reviewed in Alsos et al. 2007).All of these vectors likely could be used by fungi as well, though there is limited direct evidence concerning fungal dispersal at high latitudes (see discussion in Robinson 2001).While dispersal over long distances must have occurred in the ubiquitously distributed dominant OTUs found in our study, our data do not permit us to pinpoint the direction nor the timing of dispersal.Nevertheless, the wide distribution of EMF taxa observed in the North American Arctic and in Svalbard (Geml et al. 2012b) seem to contrast with EMF distribution patterns outside the Arctic, which often show strong phylogeographic patterns with particular taxa restricted to continents or sub-continental regions (e.g., Taylor et al. 2006, Geml et al. 2012b).
EMF species richness along latitudinal gradient.-Althougha latitudinal gradient in diversity is one of the most fundamental and striking patterns observed for many macro-organisms (Hillebrand 2004), we found no evidence of a decline with latitude in EMF species richness associated with S. arctica and D.integrifolia.This contrasts with the generally observed latitudinal species decline for vascular plants and animals in the Arctic, but our findings agree with the only other published EMF study along a latitudinal gradient, which was carried out in the European Arctic (Bjorbaekmo et al. 2010).Our study also agrees with findings for bacterial communities along latitudinal gradients in and outside the Arctic, which showed no decline in species richness (Neufeld and Mohn 2005, Fierer and Jackson 2006, Chu et al. 2010).Together, these studies suggest that species richness of prokaryotic and eukaryotic microbes may be influenced by similar factors, which appear to be distinct from those affecting macro-organisms.One confounding factor, however, is that our analysis as well as that of Bjorbaekmo et al. ( 2010) are based on the number of EMF species recorded per plant rather than the number of species per habitat or unit area.Populations of S. arctica and D. integrifolia decrease in size and are increasingly fragmented with latitude.Hence, given the well established relationships between species richness and area, it may be that EMF species richness of these two plants declines with latitude when aggregated at the level of habitat or unit area.
EMF community structure across bioclimatic subzones.-Despitethe wide distribution of the most abundant OTUs in our study, EMF community structures varied amongst bioclimatic subzones.This agrees with findings from plant communities along the same gradient (Walker et al. 2011).As observed for plant communities, EMF communities from adjacent bioclimatic subzones showed the greatest similarity.Further-  more, despite the wide distribution of the most abundant OTUs, we found a nearly complete species turnover from the coldest to the warmest subzone (Appendix: Fig. A2).This result indicates that environmental conditions, particularly temperature, contribute to shaping these communities.

Factors shaping EMF communities
Host plant identity.-Thetwo host-plants, S. arctica and D. integrifolia, which grow intermingled throughout the Arctic, shared more than a third of the observed EMF and did not harbor distinct fungal communities.The high species richness associated with S. arctica and D. integrifolia demonstrates that both hosts are broadly receptive to fungal symbionts.Our findings are in accordance with Ryberg et al. (2009Ryberg et al. ( , 2010)), who observed species-rich EMF communities on S. reticulata, S. herbacea, S. polaris and D. octopetala in a subarctic alpine tundra, with no host preference.Also, EMF communities were less host-specific in alpine tundra than in subalpine forest (Kernaghan and Harper 2001).The apparent lack of host specificity seems to be more pronounced in the arctic and alpine tundra than in boreal and temperate forests and Mediterranean woodlands, where host plants belonging to the same genera (Morris et al. 2008) or order (Ishida et al. 2007) can have an important influence on EMF fungal communities (Kernaghan et al. 2003, Ishida et al. 2007, Morris et al. 2008, Tedersoo et al. 2008, Taylor et al. 2010).The lack of a host effect on the EMF community suggests that environmental factors might be stronger drivers of this community than in some other systems.
Environmental drivers.-EMFcommunities were clearly correlated with environmental factors across regional scales and our ordinations show a strong distinction between the EMF communities from Alaska versus Canada.This distinction is likely primarily due to the glaciation history of the area and its geology.The Canadian sites were glaciated during the last Glacial Maximum and are only 10,000-16,000 years old.In contrast, the sites in Alaska, which are estimated to be 500,000-900,000 years old, were not glaciated during the Pleistocene and were part of the Beringian refugium (Raynolds and Walker 2009).Further, the sites in Canada and Alaska have a different geology resulting in different parent material and mineral contents of the soils.The sites in Canada are glacial tills and clays derived from sedimentary rocks and marine shale and are correlated with high Mg 2þ contents.In contrast, the sites in Alaska are loess, which is derived from limestone, and are correlated with higher Ca 2þ contents (reviewed in Walker et al. 2011).
EMF communities are often affected by complex site and soil factors, such as bedrock or parent material.Effects of bedrock chemistry, manifested in different pH and nutrient availability, have been observed in EMF communities associated with S. arctica on Ellesmere Island (Fujimura et al. 2008, Fujimura andEgger 2012).However, in our study, the differences in parent material affected the Mg 2þ and Ca 2þ content more than pH.All of our sites were non-acidic; the pH ranged from 6.4 to 7.8, representing nonacidic tundra.Not surprisingly, then, pH explained less of the variation of EMF community structure than in previous studies.Instead, pH becomes a key driver and strong predictor in sites with a wider range of pH values.Studies that observed such strong effect of pH on community structures include plants in the Arctic (Walker et al. 2011) and bacteria and fungi across different biomes (Erland and Taylor 2002, Fierer and Jackson 2006, Toljander et al. 2006, Chu et al. 2010, Taylor et al. 2010).
Temperature appears to be a key driver for EMF community structure along the bioclimatic gradient.This parallels findings for plant communities along the same gradient (Walker et al. 2011), as well as findings for soil fungi and bacteria along an environmental gradient in Antarctica (Yergeau et al. 2007).Experimental support for the effect of temperature on EMF community structure comes also from a longterm warming experiment in the Low Arctic tundra, where ectomycorrhizal communities of Betula nana that were dominated by Russulaceae (Russula, Lactarius) shifted to a community dominated by Cortinarius after nearly two decades of warming (Deslippe et al. 2011).While in our study Russula and Lactarius are less abundant in the High Arctic, EMF communities shifted from being Tomentella dominated in subzone B to Cortinarius dominated in subzone E. To the degree to which the bioclimatic gradient in our study is a suitable analog for climate change, this suggests that long-term climate warming of the Arctic will lead to changes in EMF communities, with an increased abundance of EMF from warmer subzones to colder subzones.The highest rates of temperature increase are observed in coastal areas of the North American Arctic (Bhatt et al. 2010) where we would expect EMF communities to change first.
Another key driver of EMF community composition is nitrogen (NO 3 À , NH 4 þ ) availability, as shown in temperate and boreal forests (e.g., Toljander et al. 2006, Kjoller et al. 2012) as well as in bacterial and fungal communities in Antarctic soils (Yergeau et al. 2007).Nitrate and NH 4 þ increased from subzone A to D and were correlated with EMF community structure, extending the patterns seen in previous studies in the Arctic.In contrast to other studies, which have suggested the importance of soil moisture on EMF abundance and community structure (e.g., Erland and Taylor 2002), soil moisture under the mesic to dry conditions in our study was not correlated with EMF community structure.
In conclusion, this is the first study of EMF communities across all five bioclimatic subzones of the North American Arctic.Our study corroborates previous findings of diverse EMF communities associated with arctic dwarf shrubs.EMF communities associated with S. arctica and D. integrifolia are dominated by a few species-rich fungal families, such as Thelephoraceae, Inocybaceae, Sebacinaceae, Cortinariaceae, suggesting that these families are particularly adapted to arctic conditions.The wide distribution of most observed EMF in and beyond the Arctic, supports emerging evidence that widespread dispersal might be a common phenomenon for fungi in the Arctic.In contrast to macro-organisms, EMF species richness does not decrease with increasing latitude and harshness of the arctic climate, at least not at the scale of individual host plants, suggesting that EMF species richness is not governed by temperature.Nevertheless, as seen in plant communities, temperature is a key factor shaping EMF community structure across the bioclimatic subzones of the Arctic.EMF communities were not affected by host identity of S. arctica and D. integrifolia, but were correlated with environmental factors across a regional scale, encompassed by a complex of glaciation history, geology, soil properties, plant productivity and climate.Our study provides important baseline data to assess climate change.Using the bioclimatic gradient as an analog for climate change, it indicates that long-term climate warming may affect EMF community structures in the Arctic by causing shifts of some EMF taxa from the Low to the High Arctic.

Fig. 2 .
Fig. 2. Species richness (number of OTUs in parentheses) within detected fungal families associated with S. arctica and D. integrifolia across the North American Arctic.

Fig. 3 .
Fig. 3. Distributions of four of the ten most abundant OTUs associated with S. arctica and D. integrifolia across the North American Arctic.Black dots indicate presence of an OTU, white dots indicate absence of an OTU in the samples studied.

Fig. 4 .
Fig. 4. Linear regression of observed (Mao Tau) and estimated (Chao2) EMF species (OTU) richness associated with Salix arctica and Dryas integrifolia along the latitudinal gradient.Species richness was rarified to three root systems.EMF species richness associated with both host plants does not significantly decline with latitude (P ¼ 0.815 (MaoTau), P ¼ 0.967 (Chao2)).We excluded ER and UP due to low sampling size ( 3 S. arctica root systems).

Fig. 5 .
Fig. 5. NMDS ordination of EMF communities associated with S. arctica and D. integrifolia for (A) the NAAT and (B) the ARCTIC datasets.The biplot diagram shows variables with r 2 ¼ 0.350 (a) and r 2 ¼ 0.100 (b).Samples were coded according to location of the study sites and host plant.

Table 1 .
Description and sampling scheme of study sites across the North American Arctic, including the Canadian Arctic Archipelago (CAA) and the North American Arctic Transect (NAAT).

Table 2 .
The ten most abundant OTUs associated with S. arctica and D. integrifolia across the North American Arctic and their occurrence on different continents, climate zones and association with different host plants as reported in GenBank.}Accessionnumbers of matching sequences from GenBank are provided in Appendix: TableA1.