Soils were sampled in close vicinity to the Zackenberg Research Station, northeast Greenland (74°28′N, 20°32′W). The climate is characterized by a mean annual air temperature of around −10°C and an annual precipitation of about 150 mm. Minimum air temperatures are below −40°C, and sub-surface temperatures (5 cm below the surface) are below −18° for about 4 months per year (Elberling and Brandt, 2003). Samples were taken from 3 different sampling sites (Table S1, Figure S1). Sites 1 and 3 were characterized by frost boils and earth hummocks and dominated by a diverse community of typical tundra vegetation (e.g., Salix arctica, Vaccinium uliginosum). Site 2 was a wet fen dominated by mosses and grasses (e.g., Eriophorum angustifolium). Two to three replicate plots that were distributed within an area of approximately 50 × 50 m were sampled at each sampling site (site 1: plots A–C, site 2: plots D–F, and site 3: plots G,H; Figure S1). For each plot, 3–5 replicate soil cores were collected close to each other (max. area: 1 × 1 m) to obtain enough material from each soil horizon. Main soil horizons were identified and pooled to obtain one replicate sample. Soil material was sampled from the active layer and the permafrost layer (PF), resulting in a total of 37 soil samples (Table S2). Soil classification follows the USDA Soil Taxonomy (Soil Survey Staff, 2010). Samples from the active layer included organic and mineral topsoil horizons (O and A), mineral subsoil horizons (B), and topsoil material that was buried into deeper soil horizons by cryoturbation (Ojj and Ajj, collectively called “J” in the following). Initial soil processing included removal of living roots prior to homogenizing the soil fraction of each sample. Samples for extracting total nucleic acids were fixed in RNAlater RNA Stabilization Reagent (Ambion Inc., Life Technologies) and kept cold until further processing. Samples for the analyses of soil properties, microbial biomass and potential enzyme activities were stored in closed polyethylene bags at 4°C until analyzed.

Soil water content was estimated by drying the soil over night at 60°C and reweighing the samples. Total organic carbon (TOC) and total nitrogen (TN) contents were determined in dried (60°C) and ground samples with an EA IRMS system (EA 1110, CE Instruments, Milan, Italy, coupled to a Finnigan MAT Delta Plus IRMS, Thermo Fisher Scientific). Microbial C and N were determined using chloroform-fumigation extraction (Kaiser et al., 2010). Fumigated and non-fumigated soil samples were extracted with 0.5 M K₂SO₄ and analyzed for extracted C and N on a TOC/TN analyzer (LiquicTOC II, Elementar, Germany). Microbial C and N were calculated as the difference in C and N concentrations between fumigated and non-fumigated samples.

Potential extracellular hydrolytic and oxidative enzyme activities involved in the degradation of organic macromolecules (namely cellulose, chitin, peptides, lignin) were measured according to Kaiser et al. (2010) using microplate fluorometric and photometric assays. The enzymes under study, their acronyms, function and substrates are listed in Table 1. One gram of sieved soil was suspended in 100 ml sodium acetate buffer (100 mM, pH 5.5) and ultra-sonicated at low-energy. Potential activities of 1,4-β-cellobiohydrolase (CBH), 1,4-β-poly-N-acetylglucosaminidase (chitotriosidase, CHT), β-N-acetylglucosaminidase (NAG), and leucine aminopeptidase (LAP) were measured fluorometrically using 4-methylumbelliferyl- (MUF) and aminomethylcoumarin- (AMC) substrates (Marx et al., 2001; Kaiser et al., 2010). 200 μL of the soil suspension and 50 μL substrate (MUF-β-D-cellobioside, MUF-N-acetyl-β-D-glucosaminide, MUF-β-D-N,N′,N″-triacetylchitotrioside, and L-leucine-7-amido-4-methyl coumarin, respectively) were pipetted into black microtiter plates in 5 analytical replicates. For each sample, a standard curve with methylumbelliferyl was used for calibration of cellobiosidase, N-acetylglucosaminidase and chitotriosidase, whereas aminomethylcoumarin was used for calibration of leucine amino-peptidase. Plates were incubated for 140 min in the dark and fluorescence was measured at 450 nm emission at an excitation of 365 nm (Tecan Infinite M200 fluorimeter). Potential phenol oxidase (POX) and peroxidase (PER) activities were measured photometrically using the L-3,4-dihydroxyphenylalanin (L-DOPA) assay (Sinsabaugh et al., 1999; Kaiser et al., 2010). Subsamples were taken from the soil suspension (see above) and mixed with a 20 mM L-DOPA solution (1:1). Samples were shaken for 10 min, centrifuged and aliquotes were pipetted into microtiter plates (6 analytical replicates per sample). Half of the wells additionally received 10 μL of a 0.3% H₂O₂ solution for measurement of peroxidase. Absorption was measured at 450 nm at the starting time point and after 20 h. Enzyme activity was calculated from the difference in absorption between the two time points.

Nucleic acid extractions were conducted according to a modified bead-beating protocol (Urich et al., 2008). Soil samples were washed with DEPC-PBS (Diethylpyrocarbonate, 0.1%; 1 × phosphate-buffered saline) to remove the RNAlater. Approximately 0.5 g of soil was added to a FastPrep™ Lysis Matrix E tube (MP Biomedicals, Solon, OH, USA). Hexadecyltrimethylammonium bromide (CTAB) extraction buffer containing 5% CTAB (in 0.7 M NaCl, 120 mM potassium phosphate, pH 8.0) and 0.5 mL phenol-chloroform-isoamylalcohol (25:24:1) was added and shaken in a FastPrep Instrument (MP Biomedicals, Solon, OH, USA) at speed 5–6 for 45 s. After bead beating, the samples were extracted with chloroform and precipitated in a PEG 6000/1.6 M NaCl solution. Pellets were washed with 70% ethanol and re-suspended in molecular biology grade water. Nucleic acids were further purified using the CleanAll DNA/RNA Clean-up and Concentration Micro Kit (Norgen Biotek Corp., Ontario, Canada). Total DNA was quantified using SybrGreen (Leininger et al., 2006).

Bacterial, archaeal and fungal SSU rRNA genes were amplified as described previously (Gittel et al., 2014). Briefly, bacterial and archaeal SSU rRNA genes were amplified with the primer set Eub338Fabc/Eub518R (Daims et al., 1999; Fierer et al., 2005) for Bacteria and Arch519F/Arch908R (Jurgens et al., 1997; Teske and Sorensen, 2007) for Archaea. Fungal SSU genes were amplified using the primers nu-SSU-0817-5′ and nu-SSU1196-3′ (Borneman and Hartin, 2000). Product specificity was confirmed by melting point analysis and amplicon size was verified with agarose gel electrophoresis. Bacterial and archaeal DNA standards consisted of a dilution series (ranging from 10⁰ to 10⁸ gene copies μL⁻¹) of a known amount of purified PCR product obtained from genomic Escherichia coli DNA and Natrialba magadii DNA by using the SSU gene-specific primers 8F/1492R and 21F/1492R, respectively (DeLong, 1992; Loy et al., 2002). Fungal DNA standards consisted of a dilution series (ranging from 10⁰ to 10⁷ gene copies μL⁻¹) of a known amount of purified PCR product obtained from genomic Fusarium oxysporum DNA by using the SSU gene-specific primers nu-SSU-0817-5′ and nu-SSU1196-3′ (Borneman and Hartin, 2000). Samples, standards and non-template controls were run in triplicates. Detection limits for the assays (i.e., lowest standard concentration that is significantly different from the non-template controls) were less than 100 gene copies for each of the genes. For statistical analyses (calculation of means and deviations per site, One-Way ANOVA), samples below the detection limit have been assigned a value of half the detection limit of the respective assay.

Nucleic acid extration and bacterial and archaeal SSU rRNA amplification were performed according to standardized protocols of the Earth Microbiome Project (www.earthmicrobiome.org/emp-standard-protocols). Amplicons were sequenced on the Illumina MiSeq platform (Caporaso et al., 2011, 2012). Sequence data along with MiMARKs compliant metadata are available from the Qiime database (http://www.microbio.me/emp/, study no. 1034).

Quality filtering of reads was applied as described previously (Caporaso et al., 2011). Reads were truncated at their first low-quality base (defined by an “A” or “B” quality score). Reads shorter than 75 bases were discarded, as were reads whose barcode did not match an expected barcode. Reads were assigned to operational taxonomic units (OTUs) using closed-reference OTU picking protocol in QIIME version 1.7.0 (Caporaso et al., 2010), with uclust (Edgar, 2010) being applied to search sequences against a subset of the Greengenes database version 13_5 (DeSantis et al., 2006) with reference sequences clustered at 97% sequence identity. Reads were assigned to OTUs based on their best hit to this database at greater than or equal to 97% sequence identity. Reads without a hit were discarded. A major advantage of this procedure is that the OTUs are defined by trusted reference sequences. It furthermore serves as a quality control filter: erroneous reads will likely be discarded as not hitting the reference data set. The primary disadvantage is that sequences that are not already known (i.e., represented in the reference data set) will be excluded. Taxonomy was assigned by accepting the Greengenes taxonomy string of the best matching Greengenes sequence. Alpha and beta diversity analyses were performed on data rarefied to 17,000 sequences per sample. Principal coordinates analysis (PCoA) was used to visualize differences between sampling sites and soil horizons using unweigthed Unifrac distances, a distance measure between communities using phylogenetic information (Lozupone and Knight, 2005). To assess differences in community structure between sites, soil horizons and replicate pits, a matrix of Bray-Curtis distances of square-root transformed relative abundances of each OTU was used for permutational ANOVA (PERMANOVA). Alpha diversity measures (richness, Shannon and inverse Simpson indices and evenness) were calculated in R 2.15.0 (R Development Core Team, 2012) using the packages vegan (Oksanen et al., 2007) and fossil (Vavrek, 2011). To test for the presence of a core prokaryotic community and distinct phylotypes in buried topsoils, we acquired the core OTUs being present in all buried topsoil samples on the non-rarified data set using the compute_core_microbiome function in Qiime.

All analyses were performed in R 2.15.0 (R Development Core Team, 2012) using the packages vegan (Oksanen et al., 2007), FactoMineR (Husson et al., 2013), and missMDA (Josse and Husson, 2013). Prior to One-Way ANOVA, Shapiro-Wilk tests and Bartlett tests were used to examine whether the conditions of normality and homogeneity of variance were met by non-transformed data sets as well as transformed data sets (log transformation: soil properties, biomass and gene abundance variables, square-root transformation: OTU abundance data). One-Way ANOVA was followed by Tukey's HSD test to determine significant differences between soil horizon (O/A, B, J, and PF) for each site. Differences were considered significant at P < 0.05.

Principle component analysis (PCA) was used to summarize data on soil properties and microbial biomass. Data were standardized (zero mean and unit standard deviation), and the number of dimensions were estimated by cross-validation and missing values were imputed using the package missMDA (Josse and Husson, 2013). Canonical correspondence analysis (CCA) was applied to relate microbial activity (potential enzyme activities) to differences in community structure (i.e., changes in relative OTU abundances). Potential enzyme activities were standardized to zero mean and unit standard deviation. Pearson's product–moment correlation (R) was used to assess linkages between individual OTUs, soil properties and enzyme activities.

Soil samples were collected from three different sites in northeast Greenland including two sites with typical tundra vegetation (site 1, n = 15 and site 3, n = 9), and a wet fen site (site 2, n = 13), and from four different soil horizons at each of these sites: organic and mineral topsoil (O and A, n = 10), mineral subsoil (B, n = 8), buried topsoil (J, n = 12), and the permafrost layer (PF, n = 7). Topsoil horizons exhibited a higher moisture, TOC and TN content, as well as higher microbial biomass (Cmic, Nmic) compared to buried topsoils, mineral subsoils and permafrost samples (Table 2). Buried topsoils were characterized by higher moisture, TOC, TN, Cmic, and Nmic than the mineral subsoils they were embedded in, although these differences were not significant in all of the sampling sites (Table 2). Moisture, TOC, and TN in topsoil and buried topsoil samples from site 2 were higher than in the corresponding samples from sites 1 and 3. PCA analysis separated O and A horizons along the first principal component and placed J horizons at an intermediate position (Figure S2). B horizons and PF samples were positioned closer to each other, thus representing higher between-sample similarity than O, A and J horizons.

On average, bacterial SSU rRNA gene abundances per gram dry soil were highest in topsoil samples and buried topsoils and were two to four orders of magnitude lower in mineral subsoils and permafrost samples (Figure 1). In contrast, average fungal SSU rRNA gene copy numbers were one to two orders of magnitude higher in topsoils than in buried topsoils. Mineral subsoils did not show any significant differences in fungal SSU rRNA gene abundances when compared to buried topsoils, while permafrost samples exhibited the lowest fungal abundances. Thus, the disproportion in bacterial and fungal abundances in topsoils and buried topsoils led to highest fungal-bacterial (FB) ratios in O horizons (FB = 3.23, Figure 1), and lowest in J horizons (FB = 0.05). Significant differences in FB ratio were found between O and A topsoil horizons as well as between O and J horizons (One-Way ANOVA, P < 0.05). Archaeal gene copies followed the bacterial distribution pattern and were highest in topsoil horizons and buried topsoils. They were two to three orders of magnitude lower in mineral subsoils and permafrost samples. The fraction of archaeal SSU rRNA genes proportionally increased with depth being lowest in O and A horizons (0.01% of the total number of prokaryotic SSU rRNA gene copies) and highest in permafrost samples (0.4%).

General observations as described above were also apparent for the individual sampling sites, including highest bacterial, archaeal and fungal abundances in topsoils and buried topsoils, and lower abundances in subsoils and permafrost samples (Figure 1). While differences in bacterial and archaeal SSU rRNA gene abundances were significant at sites 1 and 3, differences at site 2 were only minor or not significant (Figure 1). Fungal SSU rRNA abundances at site 2 were significantly higher in topsoils than in any other horizon. At site 3, fungal abundance was significantly lower in permafrost samples than in topsoil and subsoil horizons as well as buried topsoils. Differences in FB ratios between horizons were not significant (Figure 1).

Similar patterns in bacterial, archaeal and fungal SSU rRNA gene distributions were found when soil organic carbon (OC) was taken into account (Figure S3). However, differences in gene abundances per gram OC between soil horizons were less pronounced when compared to calculations on the dry weight basis (represented by lower P-values or insignificant differences, Figure S3).

Twenty-seven samples covering the four different soil horizons (O and A, B, J, and PF) were subjected to Illumina tag sequencing of the V4 SSU rRNA gene region (sequence and metadata available from http://www.microbio.me/emp/, study no. 1034). Sequencing was successful for 20 of these samples (Table 3, Table S3) and yielded a total of 1,975,888 bacterial and archaeal SSU rRNA gene sequences after extensive read-quality filtering (see Materials and Methods for details) with 86% being taxonomically classified using closed-reference OTU picking. Bacterial and archaeal sequences clustered in 10,309 OTUs, representing 157 classes (147 bacterial and 10 archaeal) within 55 phyla (including 30 candidate phyla). The dominant phyla were Proteobacteria (α, β, δ, and γ), Acidobacteria, Actinobacteria, Verrucomicrobia, Bacteroidetes, and Firmicutes accounting for ~84% of all sequences (Figure 2). In addition, Chloroflexi, Gemmatimonadetes, Planctomycetes, Cyanobacteria, and Nitrospirae were present in all soil horizons, but at lower abundances (~12% of all sequences), and 38 other rare phyla (<0.5% each) were identified. Rare phyla included archaeal taxa that represented 0.3% of all sequence reads. The majority of archaeal reads was assigned to the classes Methanomicrobia and Methanobacteria within the Euryarchaeota comprising 39.8% of all archaeal reads. For site 2, sequence reads assigned to archaeal taxa comprised 0.9% of all reads obtained from site 2 samples, whereas they only accounted for 0.05 and 0.01% at sites 1 and 3, respectively (Figure 2). This finding was in line with lower archaeal SSU rRNA gene abundances at sites 1 and 3 than at site 2 (as determined by qPCR). Furthermore, the majority of sequences retrieved from site 2 was affiliated with methanogenic taxa of the families Methanosarcinaceae, Methanobacteriaceae, and Methanosaetaceae in the phylum Euryarchaoeta, whereas members of the Thaumarchaeota and several other non-methanogenic archaea dominated at sites 1 and 3 (Figure 2).

To test whether different sites and/or soil horizons harbor distinct communities, we applied PERMANOVA and PCoA. A Bray-Curtis distance matrix of square-root transformed, relative OTU abundances was analyzed by PERMANOVA in a three-factor design to test for differences between sampling sites, soil horizons and replicate soil pits. Significant differences were returned for sampling site (Pseudo-F = 10.27, p(perm) = 0.001) and soil horizon (Pseudo-F = 6.12, p(perm) = 0.007), but not for the soil pit factor (Pseudo-F = 2.42, p(perm) = 0.097). No interactions between factors (site, horizon, pit) were apparent (P > 0.05 in all instances). PCoA based on unweighted UniFrac distances was used to visualize clusters of similar communities according to sampling site, soil horizon and soil pit (Figure 3). Most significant differences were found between sampling sites (Pseudo-F = 29.96, p(perm) = 0.001), followed by soil horizons (Pseudo-F = 22.03, p(perm) = 0.002), and soil pits (Pseudo-F = 5.93, p(perm) = 0.025). Topsoil samples from sites 1 and 3 formed a distinct cluster that also included two out of four of the buried soil samples from these sites (samples A4 and G4), but excluded topsoil and buried topsoil samples from site 2 (Figure 3). In contrast, buried topsoils from site 2 formed a distinct cluster and were thus separated from the corresponding topsoil sample as well as from samples from sites 1 and 3 (Figure 3). Subsoil and permafrost samples did not show consistent clusters according to sampling site and/or soil horizon.

Species richness and diversity were generally higher in topsoils than in buried topsoils, mineral subsoil horizons and the underlying permafrost (Table 3). However, this difference was only significant at site 2 (One-Way ANOVA, P < 0.05). Samples from buried topsoils at sites 1 and 3 were highly variable in species richness and diversity resulting in statistically insignificant differences between horizons (Table 3).

Potential hydrolytic enzyme activities (per gram dry soil) were highest in topsoil samples (O and/or A), followed by buried topsoils, mineral subsoils and permafrost samples (Figure 4). This difference was significant for CBH, CHT, LAP, and to some extent for NAG (One-Way ANOVA, P < 0.05). Potential NAG activity was similarly high in topsoil and buried topsoil samples and significantly lower in mineral subsoils and permafrost samples. Potential oxidative enzyme activities were highest in samples from topsoils (O) and buried topsoils (Figure 4). However, due to the high variability between sampling sites, differences in POX and PER activities in topsoils, mineral subsoils and permafrost samples were not significant (One-Way ANOVA, P > 0.05). Potential POX and PER activities were significantly higher in active layer samples from site 2 compared to the corresponding soil samples from sites 1 and 3 (One-Way ANOVA, P > 0.05, Figure 4, Table S4). In contrast, only few significant differences were found for the potential hydrolytic enzyme activities in samples from the same soil horizon at the different sites (Figure 4, Table S4).

When potential enzyme activities were calculated per gram organic carbon (OC), differences between horizons became smaller and less significant for hydrolytic enzymes Figure S4). For oxidative enzymes, the pattern changed to mineral subsoils (and permafrost samples) exhibiting a higher potential activity on a per gram OC basis than topsoils and buried topsoils. Differences for same soil horizons at the different sites got less pronounced for oxidative enzyme activities and did not show a clear pattern for hydrolytic enzymes (Table S5).

CCA was used to examine whether potential enzyme activities were correlated to changes in community structure (Figure 5). The first three axes of the model explained 77% of the total variation within the OTU abundance data in relation to the six enzyme activities measured. Samples from site 2 were separated from samples from sites 1 and 3 along the first axis (CC1). POX and PER activities exerted the strongest positive influence on axis 1. It was furthermore indicated that the spatial distribution of several classes within the Actinobacteria (e.g., class Actinobacteria) and the Bacteroidetes (e.g., class Bacteroidia), as well as putatively anaerobic members of the Firmicutes (Clostridia), the Chloroflexi (Anaerolinea, Dehaloccoidetes), and methanogenic Euryarchaeota (Methanomicrobia) were positively correlated with potential POX and PER activities (Figure S5).

Sequences matching putative anaerobic taxa (e.g., fermenters, sulfate reducers, metal reducers, methanogens) were highly abundant in the wet fen site accounting for up to 38% of all sequences in the O horizon (Figure 6). While the relative contribution of these taxa to the total community was similarly high in all horizons at site 2 (except for the permafrost sample), their proportion increased with depth at sites 1 and 3 (Figure 6). Here, there relative abundance increased from 2.32% in the O horizons to 9.8% in the permafrost (averaged over all samples from the respective horizons at these sites). Methanogenic archaea were particularly abundant at the wet fen site (0.5–1.3% in the active layer horizons). A high diversity of aerobic methane-oxidizing bacteria (MOB) was found at all sites and in all horizons and included MOB of type I (Methylococcaceae) and type II (Methylobacteriaceae, Methylocystaceae, Beijerinckaceae), and the recently described Methylacidiphilae within the phylum Verrucomicrobia (Table S3). Type I MOB were more abundant than type II MOB in mineral soil horizons in the active layer (38 vs. 18% of all MOB, respectively). Type II MOB predominated in topsoils and buried topsoils (38 and 46% of all MOB, respectively) together with a large fraction of putative methanotrophic Verrucomicrobia (52 and 40% of all MOB, respectively).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.