Insights into functional bacterial diversity and its effects on Alpine bog ecosystem functioning.

Plant-associated bacteria are important for the growth and health of their host, but little is known about its functional diversity and impact on ecosystem functioning. We studied bacterial nitrogen fixation and methane oxidation from indicator Sphagnum mosses in Alpine bogs to test a hypothesis that the plant microbiome contained different functional patterns depending on their functions within the ecosystem. A high abundance and diversity of nitrogenase genes were detected, mostly specific for each Sphagnum. In contrast, methanotrophs formed highly similar patterns despite a high abundance and diversity of methane monooxygenase genes. Our hypothesis was supported by these contrasting functional patterns together with the result that the Sphagnum sporophyte contained a high proportion of specific diazotrophs (45.5%) but no potential methanotrophs. While essential for plant growth under nutrient-limited conditions, nitrogen-fixing bacteria were highly specific and transferred with the sporophyte unlike the ubiquitous methanotrophs which are important for the climate-relevant ecosystem itself.

Bog ecosystems, which cover 4 Mill. km2 on Earth, have accumulated more atmospheric carbon than any other terrestrial ecosystem. Due to this extraordinary role in both carbon sequestration and net cooling effects on the global radiation balance, bogs have become essential in maintaining the global climate12. Bacterial communities living in bog ecosystems are not only involved in these important ecosystem functions, but are solely responsible for several of them3. For example, bacterial methane oxidation is crucial for the global carbon budget, and the Sphagnum-methanotroph consortium acts as a natural filter for methane456. Additionally, microbial nitrogen fixation enables plant growth under ombrotrophic, nitrogen-limited conditions in bog ecosystems7. The bryophyte genus Sphagnum, which consists of approx. 300 different species, is distributed world-wide and is the dominant component of bog vegetation8. Therefore, Sphagnum mosses have been used globally as an indicator of climate change while their microbial communities are known early indicators of ecosystem disturbances910. Recently, we have shown that living Sphagnum mosses are colonised in high abundances with specific microorganisms7111213, yet the composition, spatial distribution and transmission of the bacterial communities are well-known, less is known about the functional diversity within the microbiome.

Plant specificity of associated bacteria and fungi is a well-established phenomenon1415. Interestingly, in a study by Germida & Siciliano16, the evolutionary relationship in plant–microbe interactions was determined: old wheat cultivars were colonised by phylogenetically diverse rhizobacteria, whereas the rhizosphere of modern cultivars was dominated by fast-growing Proteobacteria. Bryophytes were the first land plants. Especially Sphagnum-dominated bog ecosystems belong to the oldest vegetation forms with more or less constant conditions for thousands of years. As a result, bog plant and microbial communities were co-exposed over long evolutionary time and developed multiple interactions. Correspondingly, an extremely high impact of the Sphagnum species was found on the structure of the microbial communities, which are transferred directly from the sporophyte (within the sporangium capsule) to the gametophyte and vice versa13. Despite this specificity, methane-oxidizing bacteria on different Sphagnum species were highly similar517. Our hypothesis for this study was that different microbial groups can show either specific or redundant pattern according to their function within the ecosystem.

We designed a new approach combining qPCR, pyrosequencing, and network analyses of functional genes to get insights into functional bacterial diversity of Sphagnum mosses and its effects on ecosystem functioning in three different Alpine bogs in Austria. While the quantitative and structural analysis of the microbiome is well-established, pyrosequencing and bioinformatic analyses of functional gene amplicon libraries are still restricted to a few, mainly methodical studies181920. We analysed key microbial genes encoding nitrogenase reductase (nifH) and particulate methane monooxygenase (pmoA) to find differences between bacterial communities of two dominant and cosmopolitan Sphagnum species, S. magellanicum and S. fallax. S. magellanicum BRID. (section Sphagnum) is typical for strongly acidic, oligotrophic and ombrotrophic habitats, whereas S. fallax H. KLINGGR. (section Cuspidata) grows in weakly acidic, more mesotrophic environments influenced by minerotrophic groundwater8. Network analysis using Cytoscape 2.8 software21 was applied to identify which functional taxa are generalists (broadly distributed across all habitats) or specialists (restricted to certain habitats but locally abundant) according to Barberán et al.22.

Results

Abundances of functional bacterial genes in Sphagnum mosses

Quantification of nifH and pmoA genes resulted in high copy numbers for both moss species (Fig. 1). In regards to nitrogen fixation, higher nifH gene copy numbers were detected for S. fallax (log10 7.0 ± 0.2 to log10 7.4 ± 0.3 copies g−1 fresh weight (FW)) in comparison to the S. magellanicum samples (log10 6.8 ± 0.1 to log10 7.0 ± 0.2 copies g−1 FW). A statistically significant difference was detected between Sphagnum spp. in sampling sites Rotmoos and Pürgschachen Moor (P = 0.004; 0.012) and in comparing copy numbers between species over all sampling sites (P = 0.0006). In regards to methane oxidation, pmoA gene copy numbers were also slightly higher for S. fallax (log10 6.3 ± 0.2 to log10 6.6 ± 0.3 copies g−1 FW) than for S. magellanicum samples (log10 6.0 ± 0.3 to log10 6.1 ± 0.2 copies g−1 FW). A statistically significant difference was shown between Sphagna from Pürgschachen Moor (P = 0.024) and testing copy numbers over all sampling sites (P = 0.002).

Figure 1

Abundance of the bacterial nifH and pmoA genes in Sphagnum gametophytes detected by qPCR.

The log abundances of nifH (a) and pmoA (b) gene copies per g fresh weight (FW) of S. magellanicum (black bars) and S. fallax (grey bars) in three sampling sites. Error bars indicate confidence intervals. Asterisks denote significant differences between moss species (*, P < 0.05; ** P < 0.01).

Deep insight into functional microbial diversity using pyrosequencing

The pyrosequencing libraries of the nifH and pmoA genes were analysed for high-throughput characterisation of diazotrophic and methanotrophic bacteria in Sphagna. The nitrogenase gene libraries were rarefied as shown in Supplementary Fig. S1. Richness estimation of the normalised datasets revealed that pyrosequencing effort attained 58.2–76.5% of the estimated richness for the clusters of 92% similarity (Table 1). The NifH clusters of 96 and 100% reflected 42.0–60.6 and 15.0–24.1% of the estimated richness, respectively. Nitrogenase diversity estimated by the Shannon diversity index (H′) was higher for S. fallax than for S. magellanicum for all sampling sites (Table 1).

Table 1

Richness estimates and diversity indices for amplicon libraries of Sphagnum samplesa

Sphagnum samplesb	Indices
Index	Clusters			Chao1			Coverage (%)			Shannon (H′)
NifH
Similarity cut-offsc	100%	96%	92%	100%	96%	92%	100%	96%	92%	100%	96%	92%
MR	618	342	122	2567	615	169	24.1	55.6	72.2	4.75	4.24	3.34
MW	702	372	145	3098	671	240	22.7	55.4	60.5	5.05	4.34	3.55
MP	636	312	108	2639	515	141	24.1	60.6	76.5	4.76	4.03	3.05
FR	724	376	178	4728	767	306	15.3	49.0	58.2	5.11	4.44	3.55
FW	778	440	216	4326	872	331	18.0	50.5	65.2	5.55	4.89	4.18
FP	846	418	194	5624	995	323	15.0	42.0	60.0	5.87	4.91	4.15

aThe number of sequences of each sample was normalised to 1658 (NifH) and 1150 (PmoA); the 16S rRNA gene amplicon library consists of 1168 sequences.

bAbbreviations specify Sphagnum species and the sampling sites: M, Sphagnum magellanicum; F, Sphagnum fallax; R, Rotmoos; W, Wasenmoos; P, Pürgschachen Moor. FSR corresponds to the sporophyte sample of S. fallax collected in Rotmoos.

cSimilarity cut-offs applied for clustering of the amino acid (NifH, PmoA) and nucleotide sequences (16S rDNA).

BLAST analysis of the NifH composition revealed an immense dominance of the amplicons related to Proteobacteria. The nifH deduced amino acid sequences were assigned to different classes and were prevailed by Alphaproteobacteria. Representative sequences of the most frequent NifH clusters showed ≥95% similarity to the NifH sequences of Bradyrhizobium, Azorhizobium, Beijerinckia, Rhodopseudomonas, Rhizobium and Methylobacterium species. Detected Betaproteobacteria were presented by the genus Burkholderia followed by Cupriavidus and Rhodoferax (≥97% similarity). NifH sequences related to Gammaproteobacteria were assigned to Stenotrophomonas and Pseudomonas genera (95–98% sequence similarity). The minor NifH clusters contained sequences 95–100% similar to Deltaproteobacteria (Geobacter spp.). A sub-dominant portion of amplicons was affiliated with cyanobacterial genera Anabaena, Tolypothrix and Aphanizomenon (≥95% similarity). According to the NCBI database, minor part of the retrieved sequences showed 98–100% similarity to the NifH sequences from methanotrophic bacteria of genera Methylocystis, Methylomonas, Methylocella and Methylosinus.

Rarefaction analysis of the particulate methane monooxygenase gene libraries resulted in similar saturation profiles for all Sphagnum samples with the exception of S. fallax from Pürgschachen Moor (Supplementary Fig. S1). While normalised datasets reached complete saturation at 82% similarity (Table 1), the PmoA clusters of 92 and 96% covered 17.4–24.7% and 46.3–94.5% of estimated richness, respectively. Except for the FP library (S. fallax, Pürgschachen Moor), methane monooxygenase diversity as estimated by the Shannon index was highly similar between Sphagnum spp. (Table 1).

PmoA diversity explored by BLAST analysis comprised type I and type II methanotrophs in an equal abundance ratio. Representative sequences of the PmoA clusters showed 89–98% sequence similarity to genus Methylomonas (type Ia), while the other group of PmoA sequences had 95–99% similarity to Methylocystis spp. (type II). Interestingly, the pmoA gene library of S. fallax from Pürgschachen Moor consisted solely of the amplicons related to the type II methanotrophs.

Comparison of the functional microbial patterns

Nitrogenase gene (nifH) amplicon libraries were compared using abundance-based Sørensen indices corrected for unseen species23. By applying 96% and 92% cut-offs, amplicon libraries formed S. fallax- and S. magellanicum-specific patterns at 55% and 67% similarity levels respectively (Fig. 2). According to the statistical analysis, 23 out of 118 examined NifH clusters (92%) showed a significant difference (P < 0.05) between Sphagnum spp. Network analysis demonstrated the distribution of the specific NifH patterns between S. fallax and S. magellanicum among all sampling sites (Fig. 3). The S. magellanicum-abundant NifH pattern consisted of sequence clusters affiliated with Methylobacterium sp., Beijerinckiaceae, Beijerinckia spp. and Azorhizobium spp. Conversely, the S. fallax-abundant NifH pattern contained clusters assigned to Burkholderia spp., Rhodoferax sp., Cupriavidus sp. and Bradyrhizobium sp.

Figure 2

Unweighted pair group method with average linkages (UPGMA) analysis of nifH and pmoA gene amplicon libraries.

UPGMA dendrograms were constructed using Chao's corrected Sørensen matrices for NifH (a–c) and PmoA (d–f) datasets of S. magellanicum (M) and S. fallax (F) from three sampling sites (R = Rotmoos, W = Wasenmoos, P = Pürgschachen Moor). Similarity cut-offs applied for clustering of datasets are indicated for each dendrogram correspondingly.

Figure 3

Profile clustering network analysis of the microbial NifH and PmoA clusters in Sphagnum samples.

NifH clusters of 92% similarity (a) and PmoA clusters of 93% similarity (b) were examined for two-fold abundance increase in Sphagnum spp. in three sampling sites. According to the results, single clusters were connected to profiles SM (abundant in S. magellanicum), SF (abundant in S. fallax) and/or COM (no change). The width of the node connections defines the number of sampling sites. The size of each node is proportional to the cluster abundance. The colour and label of the nodes specify phylogenetic and taxonomic affiliation of the clusters correspondingly. NifH network shows only clusters statistically different (P < 0.05) between Sphagnum species.

Comparison of the pmoA gene amplicon libraries revealed a high similarity between both bryophytes. With applied cut-offs, S. fallax- and S. magellanicum-originated libraries formed a common pattern of 99% similarity (Fig. 2). Statistical analysis of 23 PmoA clusters (93%) resulted in a statistically significant difference (P < 0.05) between moss species for 2 minor clusters. Network analysis revealed a general distribution of the PmoA pattern among Sphagna in all sampling sites (Fig. 3).

Linking taxonomic and functional microbial patterns of sporophyte and gametophyte

The 16S rRNA gene amplicon library of the S. fallax sporophyte generated with universal bacterial primers was analysed to reveal cross-links between taxonomic and functional diversity of the moss microbiome. Rarefaction analysis revealed complete saturation of the sporophyte library at the phyla level (Supplementary Fig. S2). At the species and genera level, a pyrosequencing survey achieved 54.8% and 65.2% of estimated richness, respectively (Table 1).

Classified 16S rRNA gene sequences were distributed among five bacterial phyla: over-dominant Proteobacteria (82.8%) followed by sub-dominant Bacteroidetes (9.4%) and Acidobacteria (7.5%) with only a few sequences affiliated to Actinobacteria and OP10. Bacteria associated with sporophyte were assigned to 37 genera (Fig. 4). According to the NCBI database, 12 detected genera (45.5% of the total diversity) comprise bacteria known for nitrogen fixation. The most dominant nifH-holding taxa were presented by Burkholderia (35.1%), Herbaspirillum (2.6%), Pseudomonas (2.0%), Sphingobacterium (1.9%), Leptothrix (1.2%) and Rhizobium (1.0%). Obtained results coincide with the network analysis: NifH clusters related to Burkholderia spp. formed a specific pattern in S. fallax gametophytes. Interestingly, no methanotrophic bacteria were detected within the sporophyte library.

Figure 4

Taxonomic classification of bacteria associated with S. fallax sporophyte.

16S rRNA gene amplicon library (FSR) was classified at genus level with confidence threshold of 50%. Only classified reads were designated (1003). Pie chart (a) represents relative abundance of the dominant genera. Genera not reaching 1% of relative abundance are assembled in the group Other and are shown in detail as a bar chart (b). Genera depicted with asterisks comprise bacteria carrying nifH gene (according to the NCBI database).

Discussion

Plant survival in nutrient-poor or ombrotrophic bog ecosystems depends on their associated bacteria. We examined two important functional groups of Sphagnum-associated bacteria in three Alpine Sphagnum-dominated bog ecosystems and found two distinct ecological patterns: diazotrophic specialists and methanotrophic generalists. In order to investigate the reasons underlying these different patterns detected by pyrosequencing of bacterial nifH and pmoA genes and corresponding network analyses, our hypothesis for this study was that different microbial groups can show specific or redundant patterns according to their function within the ecosystem. The obtained results support this hypothesis.

Nitrogen acquisition is an essential process for all forms of life, including plants24. Under nutrient-depleted conditions, bacterial nitrogen fixation could serve as the only source of N enabling plant germination and growth2526. The investigated moss species harbored highly specific nitrogen fixers which are essential for Sphagnum growth. Previous research revealed the transfer of the core microbiome from the moss sporophyte to the gametophyte13. In this study, the microbiome of the Sphagnum sporophyte contained a high proportion (45.5%) and diversity of potential nitrogen-fixing bacteria that could explain the specificity of diazotrophic communities associated with matured plants. These facts underline the importance of bacterial nitrogen fixation for Sphagnum mosses themselves.

In contrast to nitrogen fixation, methane oxidation is an important function for the whole bog ecosystem, but less essential for each member of the higher plant and moss community. Correspondingly, we detected 99% similar PmoA patterns associated with the moss gametophytes and an absence of potential methanotrophic bacteria in the sporophyte-transferred microbiome. This phenomenon can be explained as follows. Although moss-associated methanotrophs can supply the host plants with carbon dioxide4, peat mosses acquire carbon primarily from autotrophically fixed atmospheric CO2 and their carbon source is not restricted to microbial activity. On the other site, methanotrophic bacteria are able to fix atmospheric N227 and could also serve as nitrogen source for Sphagna. However, diazotrophic growth of methanotrophs is strongly sensitive to oxygen and therefore the maximum population density and activity in N-depleted bog ecosystems are located below the water table28. The investigated species S. fallax and S. magellanicum belong both to the group of non-submerged species, and therefore this overlap is unlikely. Finally, methanotrophs are dispersed by water within the bog ecosystem and thus their composition was already described as non-specific and ubiquitous517.

In our study, we applied high-throughput molecular approaches - pyrosequencing and quantitative PCR - to resolve the functional diversity of the Sphagnum microbiome. However, pyrosequencing and bioinformatic analysis of functional gene amplicon libraries are restricted to few currently published works181920 and is therefore necessary to critically assess data and conclusions. Primer selection played a crucial role for detection and comparison of diazotrophic and methanotrophic microbial patterns with pyrosequencing. A set of applied nifH gene-specific primers was shown to cover the broadest spectrum of nifH sequences29. Conversely, existing pmoA-targeting primers displayed only modest coverage, which should be considered during data analyses. Thus, the unveiled particulate methane monooxygenase diversity contained type II and type Ia methanotrophic bacteria which were ubiquitously detected through all moisture conditions, while type Ib methanotrophs specific to habitats with a high water table remained uncovered30. Interestingly, neither pmoA nor nifH gene pyrosequencing libraries contained type II methanotrophs of the genus Methylocapsa. Namely, Methylocapsa acidiphila was originally isolated from the Russian Sphagnum peat bog and exhibited acidophilic lifestyle31. When tested in silico, Methylocapsa spp. were out-selected by the pmoA gene-targeting primers and remained undetected by the pyrosequencing of the nifH gene amplicons. Furthermore, focusing on the pmoA gene diversity, we disregarded methanotrophic bacteria which use only soluble methane monooxygenase encoded by mmoX gene such as genera Methylocella and Methyloferula detected in bog environments3. In addition, a much higher diversity of methanotrophs is predicted to exist and newly described methanotrophic bacteria must be implemented in the primer design32. In spite of the partial diversity uncovered, detected Methylocystis and Methylosinus genera were shown to be dominant methanotrophs associated with living Sphagna1933.

In this study, we observed statistically greater abundances of the diazotrophic specialists as well as methanotrophic generalists in S. fallax than in S. magellanicum. This difference underlines the influence of a-biotic parameters that shape the microbiome1334. S. magellanicum lives in strongly acidic, oligotrophic, and ombrotrophic habitats, whereas S. fallax grows in weakly acidic, more mesotrophic environments influenced by minerotrophic groundwater. As a result, S. magellanicum would need more nitrogen to survive than S. fallax. Life in the ombrotrophic bog bulks is very limited, not only due to lack of nutrients also because of the extreme pH (2–4.5), which is a main driver of microbial communities35 and can explain these differences.

Recently, a high specificity of the structure of the plant-associated microbiome at cultivar level was described15. Our study shows that different functional groups are characterised by a different degree of specificity. For the first time, we found evidence for different functional pattern of plant-associated bacteria. Moreover, we found that their ecosystem function is an important driver of functional diversity.

Methods

Sampling design

To analyse differences between two Sphagnum species, S. magellanicum Brid. (section Sphagnum) and S. fallax H. Klinggr. (section Cuspidata) were selected. Both bryophytes belong to the typical and cosmopolitan vegetation in peat bogs8. Adult gametophytes of moss species were sampled in three different natural Alpine bogs in Austria in September 2009 (Supplementary Table S1). From each of the three investigated bogs, four single replicates in at least 10 m distance per Sphagnum species consisting of 15–20 plantlets were collected and stored separately. The living green parts of the plantlets were placed into sterile plastic bags and transported to the laboratory. S. fallax plants forming sporophytes were solely detected in the Rotmoos bog. Enclosed spore capsules of S. fallax were collected and processed separately. In general, S. magellanicum sporophytes are rarely found.

Isolation of total-community DNA

The microbial fraction associated with moss gametophytes was extracted as previously described13. In short, 5 g of plant material were physically disrupted with a sterile pestle and mortar and resuspended in 10 ml of 0.85% NaCl. 2 ml of suspension were then centrifuged (16500 g, 20 min, 4°C). The obtained pellet was used for isolation of the total-community DNA with the FastDNA® SPIN Kit for Soil (MP Biomedicals, Solon, OH, USA). For mechanical lysis, the cells were homogenised twice in a FastPrep® FP120 Instrument (Qbiogene, BIO101, Carlsbad, CA, USA) for 30 s at a speed of 5.0 m sec−1 and treated according to the protocol. Extraction of bacteria from the surface-sterilised sporophyte of S. fallax was carried out as previously described13. Final aliquots of the total-community DNA were further subjected to PCR-based approaches.

Quantitative real-time PCR

Quantification of the microbial nifH and pmoA genes was conducted with primer pairs nifH-F/nifH-R and A189f/A621r as referenced3637. Standards were generated by cloning the respective gene fragments from Erwinia carotovora subsp. atroseptica SCRI1043 and Methylosinus sporium ATCC35069 into the pGEM-T® Easy Vector (Promega, Madison, WI, USA). Cloned fragments were re-amplified with vector-specific primers USP (5′-GTAAAACGACAACCAGT-3′) and RSP (5′-CAGGAAACAGCTATGACC-3′) and purified using Wizard® SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA). To eliminate inhibitory effects of co-extracted substances, total-community DNA was diluted to 1:25 (data of preliminary experiments). Target genes were amplified using KAPA SYBR® FAST qPCR Kit (Kapa Biosystems, Woburn, MA, USA) in 10 μl reaction mixtures. Each measurement was performed two times in three independent runs on the Rotor-Gene 6000 (Corbett Research, Mortlake, Australia). Specificity of amplicons was confirmed by melting-curve analyses and gel-electrophoresis of the PCR products. Concentrations were calculated to copy number per g fresh weight and statistically analysed.

Deep-sequencing and bioinformatic processing

The diversity of methane-oxidizing and nitrogen-fixing bacteria associated with Sphagna was thoroughly investigated using a barcoded pyrosequencing approach. The total-community DNA of Sphagnum gametophytes was amplified with the set of nifH gene-specific primers nifH3/nifH4, nifH1/nifH2 with the nested approach38 and pmoA gene-specific primers A189f/A650r193940 using Taq-&Go™ Ready-to-use PCR Mix (MP Biomedicals, Solon, OH, USA). Total-community DNA of the S. fallax sporophyte was studied using universal bacterial primers 799f/1492r4142. Duplicate PCR products from all templates were purified with Wizard® SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA). Amplicons derived from the same Sphagnum sp. and sampling site were pooled together and subjected to the pyrosequencing using the Roche/454 GS FLX and FLX+ Titanium platforms executed by LGC Genomics (Berlin, Germany) and GATC Biotech (Konstanz, Germany) correspondingly.

Amplicon libraries of the nifH and pmoA microbial genes were studied using the FunGene Pipeline of RDP server (http://fungene.cme.msu.edu/FunGenePipeline/)43 with the parameters previously described1820. For the nifH amplicon libraries, primer sequences were trimmed and reads of low quality and shorter than 200 bp were removed. Filtered reads were translated into amino acid sequences and clipped at 60 aa. Further analyses were carried out on amino acid sequences. Rarefaction curves were generated for clusters with 100%, 96% and 92% similarity cut-offs. To calculate richness estimates and diversity indices, datasets were normalised to the same number of sequences using an in-house developed Perl script (10-times random re-sampling followed by subset formation). NifH clusters of 92% similarity (≥10 sequences) were further explored by statistical and network analyses. Representative sequences for each cluster were aligned using the NCBI algorithm TBLASTN.

The amplicon libraries of the pmoA gene were similarly processed. Trimmed reads of high quality and ≥400 bp were translated into amino acid sequences and clipped at 131 common positions. Rarefaction curves were obtained for clusters of 100%, 93% and 82% similarity. Richness estimates and diversity indices were calculated for the datasets normalised to the identical sequencing depth. The PmoA clusters of 93% similarity (≥10 sequences) were characterised through alignment of the representative sequences and subjected to further analyses.

The 16S rRNA gene amplicon library of the S. fallax sporophyte was pre-processed as previously described1213. Rarefaction analysis and richness and diversity estimation were performed using tools from RDP Pipeline (http://pyro.cme.msu.edu)43. Phylotype clusters were defined with 97%, 95% and 80% similarity cut-offs corresponding to the levels of the species, genera and phyla, respectively4445. Read classification was performed using the UCLUST pipeline integrated into the web interface SnoWMAn version 1.11 (https://epona.genome.tugraz.at/snowman/)46 with 50% confidence threshold.

Statistics

The nifH and pmoA gene copy numbers were statistically analysed using PASW Statistics 18 software (SPSS, Chicago, IL, USA). The data was tested for normal distribution by Q-Q plots and the Shapiro-Wilk test. Homogeneity of variances was checked using the Levene test. The significance of the difference between S. magellanicum and S. fallax in each sampling site and for all sampling sites was calculated using a t-test with independent samples.

To detect differentially abundant microbial clusters associated with Sphagna, NifH (92%) and PmoA (93%) clusters with ≥10 sequences were studied using Metastats web interface47. Populations of S. fallax (MR, MW, MP datasets) and S. magellanicum (FR, FW, FP datasets) were examined using a combination of the nonparametric t-test, exact Fisher's test and the false discovery rate with 1000 permutations. P-values were determined for each cluster correspondingly.

Profile clustering network analysis

Network analysis was performed to visualise functional microbial patterns and to compare their abundance between Sphagnum species across the three sampling sites48. NifH (92%) and PmoA (93%) clusters with a cumulative read change of ≥5 sequences were examined for an abundance ratio between S. magellanicum and S. fallax in each sampling site. If the ratio of values exceeded 2.0, the clusters were assigned to the specific profile (abundant in S. magellanicum or S. fallax). To refine NifH network analysis, we considered only clusters statistically different between Sphagnum species (P < 0.05). Visualization of the network was performed with the open source software Cytoscape 2.821.