Late season alpine snows are often colonized by psychrophilic snow algae that may provide a source of nutrients for microbes. Such late season snows are a harsh environment, but support a diverse and complex fungal community. We used culture independent methods (Illumina MiSeq) to test if the presence of snow algae influences fungal communities. We compared algae-colonized snows to adjacent (3 m distant) noncolonized snows in a paired experimental design. Our data indicate that several fungi are locally enriched in algae colonized snows. Although many such fungi were basidiomycetous yeasts, our analyses identified a large number of snow-borne members of phylum Chytridiomycota. While the ecology and function of these Chytridiomycetes remain unclear, we hypothesize that their enrichment in the algal patches suggests that they depend on algae for nutrition. We propose that these chytrids are important components in snow ecosystems, highlighting the underestimation of their diversity and importance. Taken together, our data strongly indicate that fungal communities are heterogeneous in snow even among adjacent samples. Further, fungal and algal communities may be influenced by similar environmental drivers resulting in their co-occurrence in snow.
Introduction
Earth's cryosphere is composed of solid water persisting for more than one month (Fountain et al., 2012) and includes all snow, permafrost, sea ice, freshwater ice, glaciers, and ice shelves. It plays important roles in global climate (Walsh et al., 2005) and terrestrial energy balance by influencing surface albedo (Prestrud, 2007). With the predicted increase in the mean global annual temperature, the cryosphere is declining and becoming a more endangered ecosystem (Derksen et al., 2012). The decrease in the annual persistence of late season snow has exceeded climate model projections in the Northern Hemisphere (Derksen and Brown, 2012). These changes likely impact local watershed dynamics as well as global biogeochemical cycles (Fountain et al., 2012). In addition to the consequences to local and global hydrological cycles, the cryosphere decline sets constraints on the distribution and diversity of organisms that depend on these unique habitats (Hoham and Duval, 2001). The transience of the late season snow coupled with its changing volume highlights the importance of understanding the endemic biodiversity residing in late season snow packs.
Snowpacks are a harsh environment characterized by low temperatures and intense ultraviolet irradiation that may act as mutagenic stressors. Yet, there is evidence for rich and diverse metabolically active microbial life in snow (Carpenter et al., 2000; Harding et al., 2011; Hell et al., 2013). These metabolically active communities contrast the belief of some that snow is but a passive recipient of aerially dispersed propagules whose metabolism is limited because of minimally available liquid water (Warren and Hudson, 2003). Yet, microbial activity or at least the presence of diverse microbial communities has been suggested for over a century (Hersey, 1913). These dichotomous viewpoints emphasize our present rudimentary understanding of microbial diversity in snow, particularly so for eukaryotic microbes. Consequently, snow can arguably be considered as a vast unexplored and undocumented ecosystem for microbial diversity (Larose et al., 2013).
Alpine and polar snows that persist through or linger late into the growing season often house snow-borne algal communities, frequently dominated by the red-pigmented algae Chlamydomonas nivalis or Chloromonas nivalis, both of whose taxonomies remain currently unresolved. These algae often produce colonies that are visible to the naked eye as a result of their characteristic red color. The red color of these algae is a result of the secondary carotenoid astaxanthin (Müller et al., 1998) and its fatty acid ester derivatives (Gorton et al., 2001) produced in large quantities during their diploid, zygotic stage that is also characterized by thickened cell walls that are resistant to harsh environmental conditions including repeated freezing temperatures (Hoham, 1980; Remias et al., 2005, 2010). During the warm season in small meltwater pools, the zygotes undergo meiosis, producing haploid offspring that are green, are metabolically active, and multiply asexually. Late in the season, when nitrogen is more limiting, these organisms develop into flagellated sexual gametes that mate producing new zygotes that can survive the next season's cold temperatures and snow (Müller et al., 1998).
Although the red C. nivalis colonies may be visually dominant, snow algal communities often consist of many species representing Chlorophyceae (Fujii et al., 2010; Remias et al., 2010). In addition, snow colonized by algae houses a broad range of fungi and bacteria that may be specifically adapted to grow in such environments (Hodson et al., 2008; Naff et al., 2013). Some evidence suggests syntrophic relationships between the snow algae and bacteria or fungi. In the most comprehensive microscopic examination to date, Weiss (1983) described the ultrastructure of the snow alga C. nivalis and repeatedly found encapsulated gram-negative bacteria on the surface of the zygotic resting stage. Weiss (1983) posited that the microscopic observations suggested syntrophy, as no similar bacteria were present in the adjacent snow without algal colonization. Similar syntrophisms have been suggested for snow algae and fungi (Kol, 1968; Hoham et al., 1993). In such syntrophic or loose symbiotic relationships, algae-associated microorganisms may utilize dissolved organic carbon (DOC) excreted by the algae. This is indirectly supported by recent studies of Chlamydomonas reinhardtii, suggesting that algae excrete large amounts of carbon (Kind et al., 2012). Algae, in turn, may benefit from the “shade” provided by the microorganisms, thus buffering the algae against the harsh environmental conditions (Light and Belcher, 1968; Hoham and Duval, 2001; Remias et al., 2005). However, Kol (1968) argued that some of these fungi might simply parasitize the algae rather than form mutually beneficial associations.
To our knowledge, studies of snow-borne fungi and bacteria are rare and limited primarily to select fungi such as “snow molds.” Generally, these “snow molds” are filamentous and thrive on organic substrates in the snow-soil interphase, but are not active in the snow itself (Robinson, 2001; Matsumoto, 2009). Recently, Naff et al. (2013) suggested that Chytridiomycota are abundant in snow and significantly influence snow food-web dynamics. Additional broad inquiries of snow-borne microbial communities indicate an abundance of microbes in snow (Harding et al., 2011) and suggest that these snow/ice-inhabiting microbes are physiologically adapted to psychrophilic environments (Gunde-Cimerman et al., 2003). Our current understanding of the general ecology of fungi in snow remains rudimentary and is based primarily on soils liberated by thawing snow (de Garcia et al., 2007) or periglacial soils (Freeman et al., 2009; Schmidt et al., 2012; Brown and Jumpponen, 2014). Studies of fungal communities associated with snow algae are usually motivated by the great abundance of fungi observed in the course of microscopic examination (Stein and Amundsen, 1967; Kol, 1968; Hoham et al., 1993).
Here we present the first high-throughput investigation of fungal communities associated with the snow algae by targeted ITS2 Illumina MiSeq sequencing. The Internal Transcribed Spacer (ITS) regions are the designated fungal barcode of life (see Schoch et al., 2012). The ITS2 region in particular has been more frequently utilized because its diminutive and largely conserved length allows for sequencing the entire region even when using more recent sequencing technologies that provide relatively short reads. We utilized a community-filtering framework (Diamond, 1975; Keddy, 1992) and explicitly test hypotheses on co-occurrence of algae and snow fungi (see Jumpponen and Egerton-Warburton, 2005). We compared paired adjacent samples with and without visible snow algae. We utilized this paired sample design because both snow algae (Müller et al., 1998) and nutrient loads (Fahnestock et al., 2000) in snow are heterogeneous and potentially confound landscape level analyses. Thus, any shifts in community-wide distributions should be detected locally (in all or most of the paired samples). Paired sample designs to interrogate community ecology are not new (Wellington, 1982; Schooley et al., 2000), but are often under-utilized in microbial community analyses (but see Taylor and Bruns, 1999; Hartmann et al., 2014). In this study, we addressed the following three questions: (1) do snow algae enrich the communities with specific fungi; (2) do snow algae shift fungal community composition; and, (3) will such shifts be consistent across larger geographic scales?
Materials and Methods
SAMPLING SITES
We sampled late season snows at six paired locations in September 2011 and August 2012 in the Glacier Peak Wilderness area, Wenatchee National Forest, Washington, U.S.A. (see Table 1 for locations and dates). Samples from each of the locations in Washington State were within 150 m of each other between 2011 and 2012. Additionally, in Colorado we sampled two paired sampling locations within Niwot Ridge Long Term Ecological Research Site ( http://culter.colorado.edu/NWT) in August 2011 and in July 2012 at two paired sampling locations within nearby Indian Peaks Wilderness area, Arapahoe and Roosevelt National Forest. Colorado sites at Niwot Ridge LTER could not be sampled in 2012 due to lack of snow in areas where 2011 samples were collected. The sampling locations were still adjacent: the maximum distance between the 2011 and 2012 Colorado sampling locations was 5.7 km. We sampled all the accessible algae colonized snows at each location, resulting in 16 paired samples (32 total). Most snow colonized by algae was inaccessible and could not be sampled (similar to Müller et al., 1998). At each sampling site, snows were collected only if the following conditions were met: (1) there were no signs of foot traffic or other anthropogenic disturbance; (2) there was an adjacent patch of uncolonized snow (based on visual assessment) within 3 m from the boundary of algae-colonized patch but within the same snow pack; and (3) there were no topographical differences between algal and nonalgal colonized snows, thus minimizing potential confounding effects of microsites. We confirmed the absence of algae in the uncolonized snows using flow cytometry.
SAMPLING PROTOCOL
Snows colonized by algae and paired noncolonized controls 3 m away from the edge of the algal colony were selected (selection criteria detailed above) and a total of five 85 cm3 volumetric surface subsamples were collected using a steel cylinder. Surface scrapings of the top 5 cm were combined into clean 1-gallon zip-top plastic bags and allowed to melt at ambient temperature. Once melted, one 100 mL sample was drawn with a sterile syringe (BD 30 mL Syringe, Dickinson and Company, Franklin Lakes, New Jersey) and passed through a 1.0 µm Nuclepore Track-Etch Membrane filter (47 mm diameter) encased in a 47 mm Swin-Lok Plastic Filter Holder (Whatman, Kent, U.K.) to collect large cells (mainly fungi and algae). The flow-through was collected into a field-sterilized container (using denatured alcohol) and passed through 0.22 µm Nuclepore Track-Etch Membrane filter (47 mm diameter) to collect bacterial cells (data not shown here). Collection and flow through containers were field sterilized with denatured alcohol between sample collections. Filter holders and syringes were used only once to minimize cross contamination between samples. After filtration, filters were stored in MoBio UltraClean Soil DNA Isolation Kit bead tubes (Carlsbad, California) with reagents S1 and IRS added to aid in sample preservation. In 2012, additional unfiltered samples were collected into sterile 15 mL falcon tubes for flow cytometric cell counts. Samples were shipped to Kansas State University and stored at -20 °C until processed. Total proportions of autofluorescent microbial constituents in snow (cells with chloroplasts relative to all cells) from the 2012 samples were estimated using flow cytometry (Guava Technologies PCA-96, Hayward, California) equipped with a 532 nm (green) excitation laser to confirm the low abundance of algae in noncolonized snows. Triplicate 10 µL samples were diluted in 250 mL 1X PBS. A combination of fluorescence emission at 675 nm (PM2, measure of chloroplast autofluorescence) and Forward Scatter (FCS) directly related to particle size were used to generate dot-plots. Using the program Flowing Software (v.2.5; www.flowingsoftware.com), boundaries of segregating autoflorescent cells along the FCS axis were manually generated based on visual clustering. The proportions of snow algae were estimated by the proportion of counts that autofluoresced but were larger than bacteria (to discriminate against cyanobacteria). Pollen grains may autofluoresce using these cytometric filters (gymnosperms have paternal ctDNA transmission and may be included) potentially overestimating the algal counts. However, pollen grain deposition in snow is generally spatially homogenous (Bourgeois et al., 2001) at the local scale and unlikely to impact the efficacy of the analyses of algal abundance. We used a paired t-test to test if the proportion of algae in the visibly colonized and uncolonized snows differed. The autofluorescent particles were on average 35 times more abundant within the algal snows than in the adjacent noncolonized snows (t = 3.25, P = 0.007). In conclusion, the algae were functionally negligible in the noncolonized snows.
TABLE 1
Sampling locations of paired algae-colonized and uncolonized snows across consecutive years. WA = Washington State, U.S.A., CO = Colorado, U.S.A. All Washington sampling locations were collected at or near the Lyman Glacier Basin.
DNA EXTRACTION AND AMPLICON GENERATION
Total genomic DNA was extracted using MoBio extraction kits according to the manufacturer's protocol with the following modifications: (1) filters were sonicated for 10 minutes (FS20; ThermoFisher Scientific, Waltham, Massachusetts) in DNA extraction tubes to dislodge any cells adhering to the filters; (2) the filter was removed and two 2.4 mm zirconia beads (BioSpec Products, Bartlesville, Oklahoma) were added into the bead tubes; and (3) particles were homogenized in a FastPrep instrument (FP120; ThermoFisher Scientific, Waltham, Massachusetts) at setting 4.0 for 60 s. The extracts were quantified (ND1000 spectrophotometer; NanoDrop Technologies, Wilmington, Delaware) and each sample was aliquoted to a 96-well plate at 2 ng µL-1 DNA concentration. PCR amplicons for Illumina MiSeq sequencing were generated using fungus specific primers to amplify the Internal Transcribed Spacer 2 (ITS2) region of the fungal rRNA gene repeat with primers fITS7 (Ihrmark et al., 2012) and ITS4 (White et al., 1990). Unique molecular identifier tags (MIDs) were incorporated to the ITS4 primer (MID-ITS4). MIDs were selected from the published Illumina MID list (Caporaso et al., 2012) and each MID-ITS4 combination was tested in silico (OligoAnalyzer 3.1; Integrated DNA Technologies, Coralville, Iowa, http://www.idtdna.com/analyzer/Applications/OligoAnalyzer) for possible hairpins and/or primer dimers at melting temperatures above 40 °C. Primers that passed this rubric were synthesized. PCR amplicons were generated in 50 µL reactions under the following conditions: 1 µM forward and reverse primers, 10 ng template DNA, 200 µM of each deoxynucleotide, 2.5 mM MgCl2, 10 µL 5× Green GoTaq Flexi Buffer (Promega, Madison, Wisconsin), 14.6 µL molecular biology grade water, and 2 U GoTaq Hot Start Polymerase (Promega, Madison, Wisconsin). PCR cycle parameters consisted of 94 °C initial denaturing step for 4 min, followed by 30 cycles of 94 °C for 1 min, 54 °C annealing for 1 min, and 72 °C extension step for 2 min, followed by a final extension step at 72 °C for 10 min. All PCR reactions were done in triplicate to control for PCR stochasticity. Negative PCR controls (sterile molecular grade water in place of DNA template) were included in each run; these controls remained free of observable amplification. Triplicate PCR products were pooled by experimental unit (total of 32) and cleaned with Agencourt AmPure cleanup kit using a SPRIplate 96-ring magnet (Beckman Coulter, Beverly, Massachusetts) following the manufacturer's protocol, except we used a 1:1 bead solution to amplicon ratio to better discriminate against small nontarget DNA. Barcoded samples were equimolarly combined so that each experimental unit was equally represented. This final pool was cleaned with Agencourt AmPure cleanup kit once more as above. Illumina MiSeq sequencing linkers were ligated onto the library and paired-end sequenced on a MiSeq Personal Sequencing System (Illumina, San Diego, California) using a MiSeq Reagent Kit v2 with 500 cycles. Ligation and sequencing were performed at the Integrated Genomics Facility at Kansas State University.
TABLE 2
The most abundant OTUs and sequence abundance for each observed phylum are provided. Nested within phylum, the most abundant order and family (genera where Incertae sedis at the family level are represented parenthetically) are provided in subsequent columns. Taxonomic representations of OTUs based on best BLASTn matches across accessioned fungi deposited in GenBank. Purported ecologies at the family (genus) level are also reported. Number of OTUs representing each taxonomic rank are provided (number of OTUs at phylum level sum to 200, and OTU counts within phylum are a subset of the total OTUs shown for phylum).
SEQUENCE ANALYSIS
Sequence data were processed using MOTHUR (v. 1.29.1; Schloss et al., 2009). The two obtained fastq (bidirectional reads) were contiged and the resultant fasta and qual files used as inputs for further sequence processing that we briefly describe below. We screened contiged sequences and required the following for inclusion: exact match to the MIDs (see Appendix Table A1 for complete list of MID sequences), at most 1 bp difference in match to both forward and reverse primers, and with a quality score of >35 over a 50 bp sliding window for the sequence. In other words, a sequence was culled if the average Q-score fell below 35 for any 50 bp window (with a 1 bp slide) or if the sequences did not match both primers or MID. Additionally, sequences were culled if they had homopolymers longer than 8 bp or contained ambiguous nucleotides. This ensured that only high-quality full-length ITS2 reads remained. ITS2 sequences were truncated to 250 bp for further analysis removing any conserved 5.8S regions, sequences shorter than 250 bp culled, and putative chimeras removed (UCHIME, Edgar et al., 2011). Remaining sequences were pair-wise aligned (Needleman-Wunsch) and the resulting distance matrix was clustered at 97% similarity using the average-neighbor algorithm. The rare OTUs (operational taxonomic units) (OTUs not among the 200 most abundant OTUs) were eliminated to focus on only the most abundant and presumably also the metabolically active taxa in these analyses. The 200 most abundant fungal OTUs represented 97.12% of all fungal sequences, therefore likely minimizing the impact of resident dormant propagules in these analyses. Randomly selected sequences representing each of the 200 most abundant OTUs were queried (BLASTn nr/nt with the exclusion of uncultured and environmental samples against GenBank) and best BLAST matches were recorded with full taxonomic string (see Appendix Table A2). Despite the use of fungal specific primers, many common OTUs were classified to nonfungal phyla: Chlorophyta (15 OTUs), Streptophyta (3 OTUs), Ciliphora (1 OTU), and supergroup Rhizaria (1 OTU). We omitted these nonfungal OTUs and appended our analyses to include only the most abundant 200 fungal OTUs. Of note is that the most abundant OTU was algal (best BLASTn match to Coenochloris sp.) and seven times more abundant than the next most abundant OTU, suggesting that primer bias was not adequate to discriminate against algal targets in samples highly enriched with phylum Chlorophyta under the reaction conditions that we used. Richness and diversity estimates were calculated iteratively—OTU richness = Sobs, Good's coverage = 1 - n l N-1 where n l is number of local singletons and N is number of sequences within each sample, complement of Simpson's Diversity = 1 - D = 1 - Σpi2 (where pi is the proportion of individuals in the ith species), and Simpson's Evenness = Ed = (1 D-1) S-1 (where S is the OTU richness at each sample and D is the Simpson's diversity index)—using 1000 iterations at a subsampling depth of 1500 sequences per experimental unit to control for biases due to unequal sampling (Gihring et al., 2011), the average of each estimator used in our analyses. This iterative approach controls for potential false positives on diversity estimates due to single subsample-based (rarefying) diversity estimates as detailed by McMurdie and Holmes (2014) by diminishing the importance of false positives as one extreme measurement has little impact on community-wide metrics after multiple iterations. Nonmetric Multidimensional Scaling (NMDS), based on a subsampled Bray-Curtis dissimilarity matrix with 1000 iterations (at 1500 subsampling depth) was used to examine fungal community composition for the first three axes (73.3% community variation, 3D stress = 0.198). To determine if the Washington sites possessed different fungal communities than the Colorado sites or if the 2011 and 2012 Washington samples differed, we used Analysis of Molecular Variance (AMOVA; PERMANOVA in Anderson, 2001). Diversity estimates and NMDS generation and analyses based on iterative subsampling were implemented in MOTHUR. Across the resolved three-dimensional NMDS space, linear (Euclidian) distances were calculated between paired algae colonized and nonalgae colonized samples and these values were analyzed using Student's t-test to test if Colorado and Washington paired samples differ in their community similarities. To test for OTU enrichment between paired algae-colonized and uncolonized snows, we used a nonparametric paired test of count data. Because of our paired design, our richness, diversity, NMDS axes scores, and OTU abundance were analyzed using a nonparametric two-tailed Wilcoxon signed-rank test (H0: Malgae=Mnon-algae,H1:Malgae≠Mnon-algae) and any significant responses were corrected for multiple comparison effects using a liberal False Discovery Rate (FDR = 0.50). Data were also analyzed using a parametric paired t-test. These analyses were consistently congruent with the nonparametric tests. As a result, we only report the nonparametric test statistics, as those rely on no assumptions on data distribution or variance homogeneity. All statistics were performed using a combination of MOTHUR and JMP (v. 10.0.2, SAS Institute, Cary, North Carolina).
The taxon assignment of OTUs to phylum Chytridiomycota (26—or 13%—of the top 200 OTUs) was challenging due to their low similarity to any vouchered or uncultured/environmental accessions. Since the low coverage and low similarity to any accessioned sequence made these BLASTn assignments to Chytridiomycota tentative, we further explored these data through a confirmatory Maximum Likelihood (PhyML) analysis. The hypervariable ITS gene regions are difficult to align globally at the phylum level and thus considered poor at discerning taxonomic relationships in phylogenetic analyses at more inclusive taxonomic ranks (order/family). Even though we used the designated barcode for identifying fungal taxa, these putative Chytridiomycota consistently failed to be assigned with satisfactory affinities. As a result, we utilized the PhyML approach as a method to confirm the phylum-level placement of the OTUs tentatively assigned within Chytridiomycota. ITS2 reads for representatives of our chytridiomycetous environmental reads (26 OTUs) were aligned with a number of accessioned and vouchered Chytridiomycota and two closely related basal phyla (Blastocladiomycota and Monoblepharidomycota) with full-length ITS2 gene regions in GenBank as well as the ascomycetous yeast Saccharomyces cerevisiae, as an outgroup. Reference sequences were selected across phylum Chytridiomycota to include all available orders of Chytridiomycota (orders Cladochytriales and Polychytriales were not included because no full-length ITS2 sequences from vouchered specimens were available in GenBank at the time of analysis). Sequences were aligned using MUSCLE (1000 iterations) as implemented in GENEIOUS (v. 5.3.4, BioMatters, Auckland, New Zealand). The alignment was trimmed manually to only include full-length ITS2 regions and a Maximum Likelihood tree was generated (100 bootstrap iterations, substitution model = HKY85) in GENEIOUS with Saccharomyces cerevisiae as an outgroup. The fasta file of the representative sequences of OTUs used for this confirmatory PhyML analysis along with reference sequences used are available as a Appendix File A1 (01-10 pp AAAR-47-4-11_A1.pdf) (available as an open access file accompanying the online version of this article).
Results
SEQUENCE DATA CHARACTERIZATION
We obtained more than 11 × 106 sequences from our MiSeq library (paired-end fastq files are deposited in Sequence Read Archive at NCBI: SRR1104197). Of these, 1,466,702 full-length ITS2 sequences were retained after quality control. After the removal of nontarget, mainly algal OTUs (363,005 sequences) and OTUs that were not among the 200 most abundant, 466,243 sequences remained with unequal sequence counts per experimental unit (range 1971–24,964). Among the total of 7835 OTUs are OTUs that are exceedingly uncommon (3949 global singletons and an additional 3686 nonsingleton rare OTUs) and most likely represent dormant organisms and aerially deposited spores or hyphal fragments that unlikely contribute to ecosystem function, or methodological artifacts of uncertain origin (Brown et al., 2015). For this reason, we limited our analysis to the 200 most abundant OTUs (see Appendix Fig. A1) that represent more than 97% of all fungal sequences in our data set (see Appendix File A2 (01-12 pp AAAR-47-4-11_A2.pdf) for OTU × Sample community matrix) ( File A2 (01-12 pp AAAR-47-4-11_A2.pdf) available as an open access file accompanying the online version of this article).
TAXONOMIC DISTRIBUTION
The fungal communities were dominated by Basidiomycota (see Table A2 for complete taxonomic assignments for the 200 most abundant OTUs and Appendix File A3 (01-19 pp AAAR-47-4-11_A3.pdf) [available as an open access file accompanying the online version of this article] for representative sequences of the 200 most abundant OTUs). Sequence and OTU counts of the most abundant orders and families with inferred ecological roles are presented in Table 2. Cosmopolitan polyphyletic yeasts in the genus Rhodotorula (46 OTUs—23% of all OTUs and 53.30% of total sequences) dominated among Basidiomycetes, as also observed in other snow and glacier fungi surveys (de Garcia et al., 2007, 2012). The four most abundant OTUs were classified as Rhodotorula and the next most abundant OTU was classified as Cryptococcus saitoi, a commonly encountered snow-dwelling yeast with uncertain ecology. OTUs assigned to Lyophyllaceae were also common (8 OTUs) and included 2 OTUs with close BLASTn affinities to accessioned sequences of Asterophora—a genus with members known to be parasites. The remaining 6 OTUs in the family were marginally similar to genus Lyophyllum, members of which are often described as soil-borne saprobic or mycorrhizal macrofungi. Chytridiomycota were surprisingly frequent (26 OTUs—13% and ∼10% of total sequences). Of these OTUs, none had a strong similarity to any accessioned sequences based on BLASTn alignments (coverage ranged from 15% to 60% with very low total BLAST scores; see Table A2). Because of the low similarity for known accessioned Chytridiomycota, we reanalyzed these data in queries that also included uncultured/ environmental sequences. These BLASTn analyses consistently failed to match closely any accessioned sequences. Despite their low similarity to sequences in the combined global genetic databases, our confirmatory Maximum Likelihood (PhyML) analysis supported placement of the 26 OTUs with 99% bootstrap support within Chytridiomycota. Additionally, our environmental OTUs were most closely related to the soil-borne chytridiomycetous order Lobulomycetales (Appendix Fig. A2). Yet, bootstrap support within the environmental Chytridiomycota clade remains low, a likely result of use of the hypervariable ITS2 region that tends to perform poorly in analyses in more inclusive taxonomic ranks (e.g., order or family). As a result, the placement of these OTUs below phylum cannot be deduced from our PhyML analysis. Yet, most of our snow Chytridiomycota OTUs form a distinct clade suggesting that these taxa may represent a monophyletic group of snow-borne Chytridiomycetes with very little ITS2 sequence similarity to anything known. Unfortunately, our phylogenetic analyses suffered from the poorly populated databases that include no vouchered reference ITS2 sequences similar to the observed environmental chytrids.
FIGURE 1.
Operational taxonomic unit (OTU) richness in the algae-colonized snows are higher than in the adjacent paired uncolonized snows (Wilcoxon Sign-Rank test). This analysis is based on 1000 iterations to generate average richness estimates (1500 sequences subsampled per iteration) based on the 200 most abundant OTUs. Uncolonized fungal OTU richness estimates are solid symbols and algae colonized richness estimates are open symbols. Dashed lines connect paired samples.
DIVERSITY ESTIMATORS
Coverage and rarefaction (see Appendix Fig. A3) estimators indicated that the fungal community was adequately captured (Good's Coverage for colonized and noncolonized snow 0.972 ± 0.003 and 0.986 ± 0.001, respectively) in our sequencing. We find that a greater proportion of the 200 most abundant OTUs occur in the algal-colonized snow than in the snow without algae. Estimated OTU richness in the algae-colonized snow (74.98 ± 11.72; mean ±1 SD) was greater than in the adjacent, noncolonized snow (64.54 ± 11.72) across both sites and years (Fig. 1; paired two-tailed Wilcoxon Sign-Rank test; W = 58, P = 0.0013). In contrast, neither diversity (complement of Simpson's Diversity, 1 - D: 0.890 ± 0.051 for the algae-colonized and 0.865 ± 0.067 for noncolonized snows; W = 30, P = 0.130) nor evenness (ED: 0.146 ± 0.053 for algae-colonized and 0.145 ± 0.073 for noncolonized snows; W = 10, P = 0.632) differed between the paired colonized and noncolonized snow samples.
COMMUNITY DIFFERENCES
Fungal communities were resolved optimally on three NMDS axes (stress = 0.198, r2 = 0.733). Our community-wide AMOVA (PERMANOVA in Anderson, 2001) analyses indicated regional (Colorado vs. Washington) but not temporal differences in the snow fungal communities. These analyses failed to distinguish fungal communities in the colonized and noncolonized snow (F 1,30 = 1.247, P = 0.227) globally. However, paired analyses of the Axis 2 ordination loading scores (representing 57.16% of community variability) indicated that fungal communities in the algae-colonized snow were distinct from those in the paired uncolonized snows (Fig. 2; W = 48, P = 0.011). In contrast, neither Axis 1 (r2 = 6.38%) nor Axis 3 (r2 = 9.85%) distinguished between the paired algal and nonalgal snows (W = 21, P = 0.298; W = 23, P = 0.252, respectively). When the paired Axis 2 values were analyzed separately for Colorado and Washington, an interesting pattern emerged. Our Colorado locations showed no discernable difference between the fungal communities of the paired samples across all three axes (see insert in Fig. 2). The Washington samples show a contrast; a shift in fungal communities in the paired algae colonized/uncolonized snows (W = 27, P = 0.034 for Axis 2; see insert in Fig. 2). Further, our AMOVA suggested no difference across years (F 1,30 = 1.135, P = 0.286) in the Washington sites but clearly distinguished the snow inhabiting fungal communities from Colorado and Washington (F 1,30 = 8.654, P < 0.001). These analyses suggest that although our sites located in Colorado and Washington are distinct, the communities remain fairly stable over time—likely as a result of local and/or regional fungal propagule inputs.
TABLE 3
Fungal operational taxonomic units that are enriched in algal colonized snow compared to paired nonalgal colonized snow based on Wilcoxon Sign-Rank test after correction for multiple comparisons (P |W| is the two-tailed P-value between colonized and paired uncolonized snow). Best BLASTn matches and putative ecologies are also reported (EcM = ectomycorrhizal). The symbol ‘‡’ represents taxa whose best BLASTn match is extremely dissimilar to any accessioned taxa (query coverage ≤ 25% and BLAST score ≤90; see Appendix Table A2) that are likely novel fungal taxa whose ecologies remain uncertain.
We also analyzed Euclidian distances between the paired samples from algae colonized and uncolonized snows to test if they differed in the magnitude of fungal community-wide shifts between Colorado and Washington. The Euclidian distances across the three resolved axes between paired algae-colonized and uncolonized snow-fungal communities was greater in the Washington snows than in Colorado (t = 2.47, P = 0.0267; Appendix Fig. A4). This difference may be the result of a constrained shift in the fungal community associated with sampling locations or dates.
In all, 13 of the 200 most abundant OTUs were more common in algae-colonized snow (see Table 3), whereas none were enriched in the uncolonized snow as determined by paired Wilcoxon Sign-Rank analyses of OTU abundance. An additional 3 OTUs differed between paired algal-colonized and uncolonized snows but were no longer significant after controlling for multiple comparisons. The colonized snow was enriched for saprobic and putatively pathogenic OTUs. Several Rhodotorula OTUs were enriched in the algal-colonized snows, suggesting the opportunistic utilization of increased organic matter associated with these snow algae. Particularly interesting are OTUs 9, 37, and 40 that—based on our PhyML analyses—represent novel chytrids whose functions remain unknown (Fig. A2). OTUs 48 and 163 were more common in algal-colonized snow but are highly dissimilar to any known accessioned fungi (see Table A2 for full BLAST scores). Best BLASTn matches identify these OTUs within the genera Tylopilus (ectomycorrhizal) and Lyophyllum (saprobic or parasitic). Given the great dissimilarity to any accessioned taxa, these OTUs most likely represent novel taxa and/or taxa that are underrepresented in the global nucleotide repositories. Thus, further and more detailed investigation is needed to better understand the fungal communities in this environment.
FIGURE 2.
Non-Metric Multi-Dimensional Scaling (NMDS) plot of snow-borne fungi in algal colonized snow (solid symbols) and uncolonized snow (open symbols). Dashed lines connect paired algal colonized and adjacent uncolonized samples. Analysis of Molecular Variance (AMOVA) indicates that Colorado (circles) and Washington (squares) fungal communities are distinct. Insert represents paired Wilcoxon Signed-Rank test across Axis 2 (57.12% of community variability) and indicates that snow algae colonized snows possess different fungal communities than paired uncolonized snows (W = 48, P = 0.011); dashed lines connect paired samples.
Discussion
We present one of the very first deep-sequencing studies of snow-borne fungi. The late season alpine snow-packs are a declining ecosystem; climate change predictions suggest that the earth's cryosphere will dramatically decline in volume (Derksen and Brown, 2012). As a result, assessment of the biodiversity in these “endangered” ecosystems is timely and critical. There is a dearth of knowledge about snow-borne life or microbial communities and their function in snow, particularly so for fungi. Our analyses suggest strongly coupled snow-borne fungi/algae co-occurrence patterns.
Our data indicate that the presence of snow algae may act as an environmental filter (Jumpponen and Egerton-Warburton, 2005) that structures the snow-borne fungal community. Alternatively, the snow algae and co-occurring fungal communities share some yet unidentified environmental variable that determines their greater frequency in those samples. This co-occurrence is best evidenced by the shift in the Axis 2 loading scores between the paired colonized and uncolonized samples (Fig. 2). Interestingly, this shift only was significant in the Washington samples in our conservative nonparametric tests. There are two primary reasons that may underlie the lack of significance in Colorado samples. (1) Our sampling in Colorado was more superficial than that in Washington and included only two algae patches that we were able to locate within the Niwot Ridge LTER or Indian Peaks Wilderness. With so few samples, we may lack the statistical power to resolve those trends. (2) Alternatively, for logistical reasons, the Colorado snows were sampled about one month before the Washington snows in both years. This temporal difference may explain this observation; perhaps the shifts in communities are not as strong in Colorado because the communities have not had as much time to diverge early in the growing season. We speculate that the establishment of greater algal abundance later in the growing season or similar environmental tolerances may facilitate a component of the extant snow fungal community. However, our Colorado sampling preceded this community shift. Our ordination analyses partially support this explanation: the Euclidian distance between the paired samples in ordination space is smaller in Colorado than in Washington (Fig. A4), suggesting that sampling later in the season may permit a further divergence among the compared snow communities. Naturally, we cannot distinguish between the temporal effects of sampling and the confounding spatial effects in this case.
Alpine snow supports diverse fungal communities dominated by basidiomycetous yeasts whose ecologies and taxonomies are poorly understood. Yeast-dominated systems have been reported on glacier surfaces (de Garcia et al., 2012) and in periglacial soils (Schmidt et al., 2012; Brown and Jumpponen, 2014). In our analyses, the most common fungal OTU was assigned to genus Rhodotorula (46 OTUs in all; 6 of which are enriched in algal snows). These yeasts are polyphyletic, understudied, and their generic delineation is historically morphological, not phylogenetic (Toome et al., 2013). Additionally, many OTUs were placed into taxa grouped as black meristematic fungi (BMF), another polyphyletic grouping based on anamorph phenotypes. Although BMFs are overall poorly understood, they are suspected to play a large role in mineral transformations and often are resilient in harsh environments (Onofri et al., 2007). These cosmopolitan BMFs most likely utilize allochthonous organic matter such as wind-blown particulate matter common on the snow surface. These results reiterate that there is a dearth of information on psychrophilic/tolerant fungi and highlights the importance for future studies into these systems.
We expected that OTUs enriched in algae-associated snows would provide the most valuable clues on the ecology of these fungi. We initially hypothesized that the algal cells and/or their nutrient-rich exudates provide substrates that potentially facilitate syntrophy or presence of opportunistic saprobes and algal pathogens. Many of the observed enriched OTUs are suspected to have saprobic/parasitic ecologies and thus support our latter hypotheses. Additionally, several taxa that were more abundant in algae-colonized snow were either macrofungi or form ectomycorrhizal associations. The presence of an active community of macrofungi in snow is unlikely but may represent snowbank fungi that fruit at the periphery of snowpacks (Cooke, 1955; Cripps, 2009). Additionally, there may be an abundance of fungal spores representing an allochthonous introduction from surrounding areas. Unfortunately, our analyses do not permit assessing whether these yeasts or other observed fungi are metabolically active in this substrate.
The abundance of Chytridiomycota in our study was surprising. This suggests that these enigmatic fungi may be more abundant and diverse than previously thought. We have just recently begun to appreciate the hidden diversity of Chytridiomycota. Freeman et al. (2009) demonstrated that Chytridiomycota dominate high altitude periglacial soils. The snow chytrids in our study were abundant and dissimilar to any sequences accessioned to the nucleotide repositories. Our confirmatory phylogenetic analyses suggested that these fungi might represent a novel clade of snow Chytridiomycota (Fig. A2). Placement of these novel chytrids from this study at levels below phylum remains uncertain but they may belong to the early divergent snow Chytridiomycetes identified as ‘Snow Clades’ from North American and European snows (Naff et al., 2013). However, this cannot be determined because our study and that of Naff et al. (2013) targeted different rRNA gene regions, making direct comparisons impossible. Naff et al. (2013) posited that these snow Chytridiomycetes parasitize snow algae because they were common in clone libraries from algal snows. However, these hypotheses were not explicitly tested, nor is parasitism the only reasonable nutritional hypothesis. The present study also differs from Naff et al. (2013) in a very important way, Naff and co-authors only collected snow that was colonized by snow algae and sequenced shallowly, whereas we utilized paired samples with and without algae. Given that Chytridiomycota also were abundant in the snows free of algae, these Chytridiomycota may be saprobes or facultative syntrophs. Nevertheless, to capture a high abundance of Chytridiomycota is striking; Chytridiomycota tend to be infrequent in most locus-targeted community sequencing studies. This may be a result of the highly divergent and difficult to amplify ITS regions of these basal fungi (Schoch et al., 2012), and analyses of environmental DNA may get overwhelmed by more easily amplified templates (Anderson and Cairney, 2004; but see Taylor et al., 2008). Thus, even the relatively high estimates of the chytrid abundance observed here might be an underestimation. Of note, solely relying on OTU abundance as a proxy for organismal abundance may be ill advised as there is no 1:1 relationship between copy numbers and organism abundance (Amend et al., 2010). Also, different fungal lineages may be differentially abundant in environmental sequences due to a myriad of factors including primer bias and differential ITS copy number (Pukkila and Skrzynia, 1993; Porter and Golding, 2012). Yet, despite the potential poor amplification of Chytridiomycota, we found that chytrids were abundant and diverse, highlighting that they both are common and likely important in snows.
It is tempting to speculate on whether these chytrids parasitize or prey on snow algae (Naff et al., 2013) as such associations are common in algae-dominated freshwater systems (Hoffman et al., 2008; Gutman et al., 2009; Rasconi et al., 2011). Snow chytrids may also act as facultative mutualists or have obligate syntrophic relationships; there is a precedence of such relationships in other aquatic systems (Picard et al., 2013). Although our data suggest the enrichment of these communities with such fungi, they do not allow specific statements about their ecology or life strategies. However, it is most likely that these novel chytrids are major players in snow-borne fungal communities. It is also clear that snow fungi are a product of establishment from local propagule pools as the snow-borne fungal communities were compositionally distinct in Washington and Colorado. In contrast, both locations had stable fungal communities over two sampling years. It is probable that snow fungi initially establish from local propagules and the presence of snow algae facilitate their growth and metabolic activity. Many small-scale biotic and/or abiotic factors that vary spatially can select for cryotolerant communities differing in their ecological and functional attributes. Further and more detailed investigations of the snow fungal metabolic activity and community dynamics are needed to better understand them in the cryosphere.
Overall, our results indicate that snow algae and snow fungi co-occur and either share similar environmental tolerances or algae may act as an environmental filter in fungal community assembly. In the latter case, this community filtering is potentially facilitated by enrichment of saprobic and pathogenic fungi that are able to utilize snow algae directly or indirectly through their exudates. Alternatively, the enrichment of specific fungal community constituents may be an outcome of facultative syntrophic associations between algae and fungi that are engaged in loose symbioses. Further in-depth studies on the life history strategies and ecology of snow-inhabiting fungi are required to shed light into these unresolved questions. Interestingly and congruently with other studies (Naff et al., 2013), our data identified potentially novel groups of Chytridiomycota, some of which were enriched in algae-colonized snows and of undetermined functions. From these studies, it is clear that we are barely scratching the surface of the nearly unexplored cryosphere. To put it simply, snow is an ecosystem that maintains unique communities that may vanish with the declining cryosphere before we have an opportunity to understand them.
Acknowledgments
The authors thank Vera Brown, Ashlyn Jumpponen, Corin White, and Kale Lothamer for assistance in sampling. We thank anonymous reviewers of previous versions of this contribution for providing useful suggestions that improved this manuscript. We are indebted to those anonymous strangers who assisted Brown during a field injury. We are extremely appreciative to the Douglas County, Washington, Sheriff's Office for not confiscating our samples when dispatched about suspicious people in a park with gloves and syringes. Brown was supported by a National Science Foundation (NSF) GK-12 fellowship (DGE-0841414 to C. Ferguson) and the Kansas Academy of Science.
References Cited
Appendices
APPENDIX
TABLE A1
Primer sequences (upper) and multiplex Molecular Identifier sequence tags (lower) used in this study.
TABLE A2
Full taxonomic descriptions of the 200 most abundant fungal OTUs including sequence count, best BLASTn match, Max Score, Total Score, Query Coverage, E-Value, Max Identity and Accession numbers for the closet match.
Continued
Continued
Continued
Continued
Continued
FIGURE A1.
Ranked operational taxonomic unit (OTU) abundance distribution plot of all OTUs. The dashed line represents our cutoff of 200 OTUs for analyses and represents greater than 97% of all fungal sequences. Insert represents the 30 most abundant OTUs (notice the change in scale) demonstrating that few OTUs make up most of the fungal community members.
FIGURE A2.
Phylogenetic analysis (Maximum Likelihood) of putative chytrid OTUs with vouchered representative ITS2 sequences within phyla Chytridiomycota, Blastocladiomycota, and Monoblepharidomycota indicates that observed novel OTUs are nested within phylum Chytridiomycota with 99% bootstrap support.
FIGURE A3.
Rarefaction analysis of observed OTUs for algal colonized and non-algal colonized snow fungi indicate that at the 1500 sequence subsampling point (dashed line), the majority of community members have been observed as this subsample value is well past the inflection point of the curves.
FIGURE A4.
Distances of paired algal colonized and non-algal colonized axes loading score across three-dimensional space indicate that Colorado (CO) fungal communities are more similar between paired samples than Washington (WA) paired fungal communities (Euclidian distance between paired samples based on t-test).
