The biogeochemical variability of Arctic thermokarst ponds is reflected by stochastic and niche-driven microbial community assembly processes.

Shallow thermokarst ponds are a conspicuous landscape element of the Arctic Siberian tundra with high biogeochemical variability. Little is known about how microbes from the regional species pool assemble into local pond communities and how the resulting patterns affect functional properties such as dissolved organic carbon (DOC) remineralization and greenhouse gas (GHG) turnover. We analysed the pelagic microbiomes of 20 ponds in north-eastern Siberia in the context of their physico-chemical properties. Ponds were categorized as polygonal or trough according to their geomorphological origin. The diversity of bacteria and eukaryotic microbes was assessed by ribosomal gene tag sequencing. Null model analysis revealed an important role of stochastic assembly processes within ponds of identical origin, in particular for genotypes only occurring in few systems. Nevertheless, the two pond types clearly represented distinct niches for both the bacterial and eukaryotic microbial communities. Carbon dioxide concentration, indicative of heterotrophic microbial processes, varied greatly, especially in the trough ponds. Methane concentrations were lower in polygonal ponds and were correlated with the estimated abundance of methanotrophs. Thus, the overall functional variability of Arctic ponds reflects the stochastic assembly of their microbial communities. Distinct functional subcommunities can, nevertheless, be related to GHG concentrations.

Introduction

Lakes and ponds are common landscape elements of the arctic permafrost lowlands. The most prevalent water bodies of the northern tundra are small shallow ponds with a surface area under 1 km2 and a depth lower than 2 m (Verpoorter et al., 2014; Muster et al., 2017). These small ponds are disappearing in some regions like the Alaskan North slope (Riordan et al., 2006; Andresen and Lougheed, 2015), while expanding in size and numbers in others, such as subarctic Canada (Payette, 2004) and north‐eastern Siberia (Walter et al., 2006). Although the hydrology of these ponds seems to follow natural cycles of formation, expansion and drainage (van Huissteden et al., 2011), climate warming is likely to trigger an acceleration of these cycles in the Arctic regions (Vincent et al., 2013).

Given their shallowness and small size, arctic ponds are characterized by a large interface between the water and the sediments, the surrounding terrestrial ecosystems and the atmosphere. Although ice‐covered or frozen solid during most of the year, they are often nutrient‐rich and support high primary productivity (Przytulska et al., 2016) as well as active bacterial communities (Breton et al., 2009), highlighting their role as biogeochemical hotspots (Laurion et al., 2010; Deshpande et al., 2016). Arctic terrestrial ecosystems are particularly sensitive to climate change, resulting in the rapid thawing of the permafrost, thereby mobilizing the important stock (Tarnocai et al., 2009) of ancient and labile contemporary organic carbon (OC) into the shallow thermokarst ponds (Schuur et al., 2015). Although periodically acting as CO2 sinks (Laurion et al., 2010), thermokarst ponds are often CO2 and CH4 oversaturated and thus generally recognized as a significant source of greenhouse gases (GHGs) (Sobek et al., 2005; Breton et al., 2009). CH4 emissions from such ponds account for two‐third of the total emissions in regions above the latitude 50°N (Wik et al., 2016) and are strongly stimulated by atmospheric warming (Negandhi et al., 2014; Negandhi et al., 2016; Bartosiewicz et al., 2019). However, shallow ponds are still mostly neglected in estimates of regional and global C budgets (Muster et al., 2012).

Ponds in the Arctic originate from various geomorphological processes (Pienitz et al., 2008) that can lead to diverse biogeochemical properties (Breton et al., 2009). For instance, polygonal ponds form during the freeze–thaw cycles of ice wedges (Lachenbruch, 1962) and display transparent waters with benthic microbial mats (Bonilla et al., 2005). Other ponds form in depressions created by the underlying thawing permafrost and display turbid waters containing high amounts of dissolved OC (DOC) (Rautio et al., 2011). Shallow ponds also show remarkable spatio‐temporal variability in their GHG emissions (i.e. CH4 emissions ranging from 0.03 to 5.62 mmol m−2 day−1; Laurion et al., 2010). Although DOC concentration is an important factor controlling GHG emissions (Breton et al., 2009), the reasons for such disparities remain poorly understood.

The biogeochemical cycle of C in aquatic systems is mainly driven by microbial activity (Del Giorgio and Cole, 1998). The rates of C remineralization and storage depend on the interplay of abiotic conditions, such as pH and temperature (Muscarella et al., 2019). At the same time, it is unclear to which extent the microbial community composition and process rates actually reflect the local environmental conditions (Langenheder et al., 2006; Graham et al., 2016). In view of stochastic community assembly processes such as dispersal and ecological drift, the specific traits of individual microbial taxa may not be decisive for their respective occurrence patterns (Stegen et al., 2013; Langenheder et al., 2017). The study of β‐diversity allows to assess the link between the environment, biodiversity and ecosystem functioning. In this context, null model approaches are powerful tools to assess the importance of various community assembly processes. The Raup‐Crick (RC) metric represents a significance test that allows to disentangle changes in β‐diversity from variation in α‐diversity (Chase et al., 2011) while the normalized stochasticity ratio (NST) is an index that has been developed to quantify the relative importance of deterministic and stochastic community assembly processes (Ning et al., 2019).

While microbial community composition across permafrost thawing gradients in the Canadian Arctic appears to be mainly shaped by environmental filtering (Crevecoeur et al., 2015), stochastic community assembly processes can be dominant in ponds formed on organic permafrost (Comte et al., 2016). So far, only few studies have sought to integrate the variability of microbial community structure with that of physico‐chemical properties to assess the role of small thermokarst ponds in the C cycle (Vonk et al., 2015).

We investigated the bacterial and eukaryotic pelagic microbial communities in the context of physico‐chemical variability in two distinct types of shallow Arctic ponds with the following objectives: assess (i) if the microbial community structure in these ponds was primarily shaped by neutral assembly processes or determined by differences in habitat structure and/or pond origin and (ii) if the microbial assembly processes and resulting structure were reflected in the biogeochemistry of these shallow ponds in the context of GHG saturation and potential emissions.

Results

Morphometric properties and surrounding vegetation

All ponds were small with a surface area not exceeding 150 m2 (Fig. 1, Fig. S1). The polygonal ponds (site A) were on average larger than the trough ponds (site B) (Table 1) with large variability between the ponds. All ponds were very shallow, their depth ranging between 20 and 30 cm (Table 1). The active layer at the bottom of the ponds ranged from 24 to 40 cm, while the active layer in the neighbouring soil was between 17 and 30 cm. The two types of ponds were surrounded by different plant species (Fig. S2). The dominant vegetation around the polygonal ponds was Sphagnum fuscum and Carex aquatilis, while the trough ponds were surrounded by Ledum palustre, Eriophorum vaginatum, Betula nana and the moss Dicranum. The trough ponds were characterized by the presence of submerged Betula nana (Fig. 1D).

Fig. 1

a. The Kytalyk region in North‐Eastern Siberia, Russia. b. Location of the two sampling sites, situated 1.5 km apart. The star represents the camp. c. Example of a polygonal pond at site A. d. Example of a trough pond with flooded vegetation at site B.

Table 1

Comparison of morphometric characteristics and physico‐chemical properties of the 10 polygonal (site A) and 10 trough (site B) ponds (means ± standard deviation).

Property	Polygonal (A)	Trough (B)	p‐value
Elevation (m)	17.7 ± 0.92	13.5 ± 0.69	0.0007
Surface area (m2)	144.3 ± 100.4	57.9 ± 55.7	0.03
Depth (cm)	22 ± 2.5	23 ± 4.1	n.s.
Active layer depth (cm)	34.3 ± 1.9	33.8 ± 4.9	n.s
pH	5 ± 0.1	5 ± 0	n.s
Temperature (°C) a.	15.0 ± 1.1	15.9 ± 1.2	n.s
DOC (mg l−1)	21.54 ± 5.13	62.02 ± 18.44	0.0002
DN (mg l−1)	0.01 ± 0.02	0.95 ± 0.41	0.0007
DIC (mg l−1) a.	0.62 ± 0.39	0.23 ± 0.20	n.s
cDOM a320 (m−1)	41.06 ± 16.6	170.97 ± 58.5	0.0003
SUVA254 (l mg−1 m−1) a.	1.98 ± 0.40	2.81 ± 0.38	0.0009
Slope ratio	3.96 ± 0.15	3.68 ± 0.07	0.0006
POC (mg l−1)	2.45 ± 1.29	5.82 ± 4.16	0.02
δ 13C POC (‰)	−32.3 ± 2.05	−32.6 ± 1.95	n.s
dissolved CO2 (μM)	11.61 ± 10.24	26.62 ± 25.10	n.s
CO2 saturation (%)	65.7 ± 58	150.6 ± 142	n.s.
dissolved CH4 (μM)	0.76 ± 0.61	5.99 ± 7.78	0.006

a. Comparisons of these mean values were performed with 2‐sidedt‐tests. The other parameters were compared by nonparametric Mann–Whitney U‐tests.

Significant differences are highlighted in bold. DOC: dissolved organic carbon, DN: dissolved nitrogen, DIC: dissolved inorganic carbon, cDOM a320: absorption coefficient by chromophoric dissolved organic matter at 320 nm, SUVA254: specific UV absorbance at 254 nm, POC: particulate organic carbon.

Physico‐chemical characteristics

Chlorophyll‐a (Chl‐a) concentrations were on average 1.1 μg l−1 (range 0.2 to 5.2 μg l−1) across all ponds (excluding pond B05 that had an exceptionally high Chl‐a concentration of 38.1 μg l−1), without differences between the two types (Mann–Whitney, W = 54, p = 0.5). DOC and cDOM were three and four times higher, respectively, in the trough ponds than in the polygonal ponds (Table 1). The trough ponds were also enriched in POC and DN compared to the polygonal ones, while the latter contained slightly more DIC. The DOC of the trough ponds consisted of more aromatic (higher SUVA254) and higher molecular weight compounds (lower slope ratio) than that of the polygonal ponds. The δ13C isotopic signal of POC averaged −32‰ in both pond types. Differences in carbon content were also visible by colour: the trough ponds were dark brown, while the polygonal ponds were transparent. Ponds from both sites displayed pronounced variability in their water isotopic signatures (Fig. S3), and their local evaporation line deviated from the global meteoric water line. The water of all ponds was enriched in heavy water isotopes, with polygonal ponds being significantly more enriched in deuterium (t‐test, t = 3.8, df = 18, p = 0.004) but not in 18O (t‐test, t = 1.6, df = 18, p = 0.15). The ponds were acidic with a pH of five and the water temperature was uniform over the sampling period. The amount of dissolved CO2 in the ponds varied greatly and oversaturation (CO2 equilibrium at 15.5°C: 17.7 μM) was found in five trough ponds and one polygonal pond, while the others were CO2 depleted (Table 1). All ponds were oversaturated in dissolved CH4 compared to atmospheric concentrations (CH4 equilibrium at 15.5°C: 3.1 × 10−3 μM). Dissolved methane concentrations were highly variable, ranging from 0.18 μM to 27.5 μM, and they were significantly higher in the trough ponds (Table 1). The δ13C‐CH4 signatures varied among the polygonal (mean −48.9‰, range −58.4‰ to −30.3‰) and the trough ponds (mean −49.0‰, range −69.8‰ to −16.5‰). The range of variability was somewhat larger in the trough ponds but this difference was not statistically significant (F test, F = 0.28, p = 0.11). The water isotopic signal (δ2H) was correlated with the CH4 concentration (Spearman's ρ = −0.51, p = 0.04) as well as with the δ13C‐CH4 signatures (Pearson's r = 0.48, p = 0.04).

Microbial community composition

Bacterial abundances ranged from 3.9 to 10.3 × 106 cells ml−1 in the polygonal ponds and from 4.8 to 29.1 × 106 cells ml−1 in the trough ponds, with no significant difference between pond types (median values for polygonal ponds: 6.1 × 106 cells ml−1, trough ponds: 7.4 × 106 cells ml−1, Mann–Whitney, W = 34, p = 0.34). Bacterial OTU numbers in the normalized data set varied from 60 to 420 OTUs (Fig. 2A), without difference between pond types (t‐test, t = −0.46, df = 38, p = 0.64). Bacterial abundances and OTU numbers did not correlate with DOC (respectively: Spearman's ρ = 0.23, p = 0.3; Spearman's ρ = 0.13, p = 0.5). OTU numbers of eukaryotic microbes ranged from 6 to 255 OTUs (Fig. 2B), with polygonal ponds being more diverse than the trough ponds (t‐test, t = 2.57, df = 38, p = 0.01). Pond B05 displayed the lowest diversity of both, bacteria and microbial eukaryotes.

Fig. 2

Number of OTUs in ponds from the two sites: (a) bacterial communities, (b) microbial eukaryotic communities. Asterisks indicate significant differences at p < 0.001.

Bacterial communities in both pond types were dominated by Proteobacteria, Actinobacteria and Bacteroidetes (Fig. 3A). The most ubiquitous and numerous OTUs were a sequence affiliated to the genus Polynucleobacter and an OTU that matched 100% to database entries assigned to both, Variovorax (accession number: MT102304.1) and Rhodoferax (accession number: MH699237.1) in BLAST. We also noted the presence of OTUs affiliated to the Candidate Phyla Radiation. Similarity percentage analysis (SIMPER) indicated that OTUs affiliated with Variovorax/Rhodoferax, candidatus Planktoluna, Limnohabitans and Rhodoblastus preferably inhabited the trough ponds while an uncultured Burkholderiaceae (with a similar sequence to Alcaligenes), Rhodoluna, candidatus Planktophila (acl‐A) and Methylorosula were more common in the polygonal ponds (Table S1).

Fig. 3

OTU composition within different phylogenetic lineages of the (a) bacterial and (b) microbial eukaryotic communities of polygonal ponds and trough ponds. OTUs in the polygonal ponds (left side of panels) were sorted according to their contribution to all OTUs in this pond type (y‐axis), and plotted against their corresponding contribution in the trough ponds (right side of panels). The numbers in the right half of the panels report the relative contribution of individual phylogenetic lineages to the average dissimilarity between polygonal and trough ponds communities (SIMPER analysis of Bray‐Curtis distances, significantly contributing OTUs only).

Methanotrophs belonged to the genera Methylocapsa, Methylocystis, Methylosinus, Methylocella, Methylobacterium, Methylococcus, Crenothrix, Methylomonas, Methyloparacoccus, Methyloglobulus and candidatus Methylospira (Fig. S4). On average, 1.2% (range 0.02–9.1%) of the bacterial reads in the polygonal ponds and 0.2% (range 0–0.9%) in the trough pond were affiliated with these taxa, with some indication of differences in their community composition between pond types (Fig. S4). More methanotrophs (estimated abundance) were present in polygonal ponds (Mann–Whitney, W = 334, p < 0.001). Their estimated abundance was moreover negatively correlated with dissolved methane concentrations (Spearman's ρ = −0.63, p = 0.01) (Fig. 4) and positively correlated with the δ13C‐CH4 signatures (Spearman's ρ = 0.51, p = 0.04).

Fig. 4

Relationship between the estimated abundance of methanotrophs and dissolved CH4 concentration. Data were normalized by a Box‐Cox transformation (λ = −0.14). Pond B05 was excluded from the analysis due to its lack of methanotroph reads (see text).

The eukaryotic microbial communities in both pond types mainly consisted of members of the SAR (Stramenopiles, Alveolata, Rhizaria) lineage, Archaeplastida, Cryptophyceae and fungi (Fig. 3B). The difference between the communities of the polygonal and trough ponds (SIMPER) was mainly due to members of Archaeoplastida (genera Monomastix, Chaetosphaeridium) and Stramenopiles (Uroglena) that were more abundant in the polygonal ponds (Table S2). The trough ponds had higher proportions of reads from a single OTU that was 100% identical to Cryptomonas marssonii and that accounted for 21% of the separation between pond types. Pond B05 (which also exhibited unusual high Chl‐a concentrations, see above) was dominated by a genotype with 98.4% sequence identity with a Chlamydomonas proboscigera (GU117581) isolated from Svalbard. Fungi were also an important component of the eukaryotic communities, mainly represented by the class Chytridiomycetes.

β‐diversity of microbial communities across ponds

Within each pond type, the average pairwise Bray‐Curtis (BC) similarity of bacterial communities was lower than 50% (polygonal ponds: mean 33%, range 5–78%; trough ponds: mean 34%, range 3–80%). Both types of ponds showed similar variability in bacterial community structure (homogeneity of multivariate dispersions, N.perm = 9999, F = 0.007, p = 0.92). However, polygonal and trough ponds were composed of distinct communities (PERMANOVA, N.perm = 9999, F = 6.9, R 2 = 0.16, p < 0.001) (Fig. 5A), with low pairwise similarity between ponds of different types (mean 25%, range 3–77%).

Fig. 5

Non‐metric Multidimensional Scaling (nMDS) of (A) the bacterial communities (stress = 0.15) and (B) the eukaryotic microbial communities (stress = 0.14) based on Bray‐Curtis β‐diversity. The ellipses represent the 95% confidence interval of the position of each group.

The BC similarity of the eukaryotic microbial communities was low within polygonal ponds (mean 19%, range 2–51%), whereas trough ponds were on average more similar to each other (mean 71%, range 0–90%). The variability of eukaryotic community composition was similar for both pond types (homogeneity of multivariate dispersions, N.perm = 9999, F = 0.88, p = 0.36). Nevertheless, the microbial eukaryotic communities significantly differed between pond types (PERMANOVA, N.perm = 9999, F = 9.6, R 2 = 0.20, p < 0.001) (Fig. 5B), with very low similarity between pond pairs from different types (mean 0.07%, range 0–37%).

Microbial community assembly processes

Stochastic processes dominated the assembly of bacterial communities within both pond types, as indicated by NSTs > 50% and more than 80% of RC indices ranging between −0.95 and 0.95 (Table 2). The contribution of stochastic processes was similar in both pond types (bootstrap test on NST, p = 0.39). However, two lines of evidence suggested that the two sets of ponds represented distinct niches: their between‐type NST was < 50%, and approximately 30% of the between‐type RC comparisons were less similar than expected by chance. There was no clustering nor over‐dispersion of phylogenetic diversity either within or between the two sites [−2 < Beta Nearest Taxon Index (βNTI) < 2].

Table 2

Community assembly processes of polygonal ponds (site A) and trough ponds (site B) as assessed by null model based indices.

Bacteria	Polygonal (A)	Trough (B)	A versus B
NST	63.4	65.4	43.9
RC < −0.95	2.3	2.3	2
−0.95 < RC < 0.95	84.4	97.7	63
RC > 0.95	13.3	0	35
βNTI	−1.44	−0.82	−1.10

βNTI: Beta Nearest Taxon Index, NST: normalized stochasticity ratio (%), RC: modified Raup‐Crick (% of pairwise comparisons within specified range). Polygonal, trough: within sites; A versus B: between sites.

The community assembly processes of the eukaryotic microbial community appeared to be more related to the pond type. Specifically, stochastic processes dominated in the polygonal ponds and deterministic processes in the trough ponds (Table 2). However, this difference was not statistically significant (bootstrap on NST, p = 0.07). Interestingly, 35% of polygonal pond pairs were less similar than expected by chance (RC > 0.95), whereas the opposite appeared in the trough ponds, with 13% of the pairs being more similar than expected by chance (RC < −0.95). Pond type determined the microbial eukaryotic communities even stronger than the bacterial communities (NSTA vs. B = 30%), and communities originating from different pond types were less similar than expected by chance in 62% of the RC comparisons. The phylogenetic composition of the eukaryotic microbial community was neither clustered nor over‐dispersed within and between the two types of ponds (−2 < βNTI < 2). The NST indices of bacterial communities from both pond types steeply decreased if only the more widely distributed OTUs were considered (Fig. 6), and it rose again for the small subset of the most ubiquitous OTUs that were present in 18 ponds or more. No such trend was observed for eukaryotic microbes, and no eukaryotic OTU occurred in more than 16 of the ponds.

Fig. 6

Normalized stochasticity ratio (NST) of subsets of OTUs occurring in increasing numbers of ponds (n max = 20 and 16 ponds for bacteria and eukaryotes respectively).

The influence of deterministic processes on bacterial community composition was also indicated by a weak but significant correlation with the measured environmental variables (Mantel test, 9999 permutations, Spearman's ρ = 0.20, p = 0.04). The distance‐based redundancy analysis (db‐RDA) (Fig. S5A) was significant (df = 3, R 2 = 0.35, p = 0.001) with DOC, DIC and δ13C POC as the main explanatory variables and explained nearly 50% of the total inertia (CAP 1: 36%, CAP 2: 12%). The microbial eukaryotic communities were also significantly correlated with the environmental variables (Mantel test, 9999 permutations, Spearman's ρ = 0.51, p < 0.001). The db‐RDA (Fig. S5B) was also significant (df = 7, R 2 = 0.28, p = 0.001) with DOC, DIC and POC as the main explanatory variables. Only the CAP 1 explained inertia (21%).

Discussion

Pond types differ in habitat properties

The shallow ponds of the Kytalyk region, although closely located to each other, were remarkably variable in their biogeochemical properties (Table 1), as has also been observed in thermokarst ponds in other circumpolar regions (Breton et al., 2009; Shirokova et al., 2009; Rautio et al., 2011). The two pond types were set apart particularly by the quantity and optical properties of OC (DOC and POC). The OC in thermokarst ponds mainly originates from terrestrial input (Crump et al., 2003; Wauthy et al., 2018). It has a great influence on pond colour and turbidity, which in turn affects primary production and subsequent trophic processes (Fasching et al., 2014). The observed difference in DOC properties between pond types could partly be due to features of the catchment, such as the type of the adjacent vegetation or of soil, both known to influence the quantity and biodegradability of allochthonous DOC (Blom et al., 2000; Larouche et al., 2015). Indeed, the surrounding vegetation clearly differed between the sites (Fig. S2). Since the ponds were situated within the beds of two ancient thaw lakes of different drainage age, it is conceivable that the composition of the active layer and permafrost in carbon and ice might differ as well. The geomorphological origin and maturation of thermokarst ponds also influence the amount of C discharged into ponds (Pokrovsky et al., 2011; Peura et al., 2020). The banks of trough ponds appeared to be less stable than of the polygonal ones, as indicated by the presence of flooded Betula nana.

Photolysis in Arctic aquatic systems can be an important process in modifying DOC optical properties and bioavailability (Vähätalo et al., 2003; Laurion and Mladenov, 2013), especially in summer where surface waters are exposed to a 24 h light regime. The low amount of cDOM in the polygonal ponds in combination with higher DIC concentrations suggests that DOC underwent more intense photochemical degradation in these systems than in the trough ponds. However, bacterial abundances and α‐diversity were similar in both types of ponds and did not correlate with DOC quantity and type. This stands in contrast to observations from subarctic ponds (Roiha et al., 2012) and boreal lakes (Logue et al., 2012) and suggests that bacterial abundance and α‐diversity were not primarily controlled by DOC availability.

The bacterial and eukaryotic microbial communities of polygonal and trough ponds clearly differed in composition (Fig. 5), and null model analyses suggested that these differences were predominantly driven by deterministic processes (Table 2). Since the pond types were distinguished by geographic location and habitat characteristics (Fig. 1, Table 1), two scenarios could explain this observation. For one, each type of ponds might have a different local species pool as a result of dispersal limitation, which may be relevant even at small spatial scales (< 1 km) (Martiny et al., 2011). On the other hand, some studies have argued that dispersal limitation is only relevant at larger spatial scales and in physically well‐separated landscapes (Green et al., 2004; Comte et al., 2016). Assuming a minor role of dispersal limitation, the observed occurrence pattern of taxa in different pond types would be best explained by niche‐driven processes. In the case of bacteria, this conclusion is also supported by db‐RDA analysis (Fig. S5), with DOC explaining an important part of the variance in community composition. However, an NST index close to 50% and the 63% of the RC indices consistent with the null model both indicate that stochastic processes still played an important role in shaping bacterial communities between pond types. Deterministic factors between pond types were more important in the micro‐eukaryotic assemblages (Table 2). However, the link with the measured environmental parameters was weaker (Fig. S5), indicating that the relevant niche parameters were not accounted for.

Stochastic processes dominated the assembly of bacterial communities within both pond types and of eukaryotic microbes in polygonal ponds. Priority effects (e.g. impact of colonization order) and ecological drift (e.g. random changes in relative abundances of species) are major stochastic processes (Zhou and Ning, 2017) that might play an important role in the community composition of such ponds. Ponds represent small habitat patches and are seasonally ephemeral (they freeze solid during winter). This speaks for frequent pond‐specific colonization events, favouring priority effects (Fukami, 2015).

In the trough ponds, the NST index of eukaryotic microbial communities pointed to a dominance of deterministic assembly processes, and 15% of the RC indices in these ponds were more similar than expected by chance (Table 2). Thus, eukaryotic microbes might have experienced species sorting by pond characteristics that were not explicitly determined in this study, such as limited light availability for phototrophs (Thrane et al., 2014) or a stronger top‐down‐control by daphnids (Zöllner et al., 2003; Sarnelle, 2005) that were only present in trough ponds (A.L.M., personal observation).

Community assembly processes may vary for rare and abundant taxa (Liu et al., 2015; Nino‐Garcia et al., 2016). We observed a clear decrease of the NST index in subcommunities of those bacterial OTUs that occurred in an increasing number of ponds (Fig. 6A). This suggests that the stochastic overall assembly pattern of the bacterial communities was primarily shaped by taxa that were only present in a single or a few systems. Bacteria are inoculated into the ponds mainly from the surrounding soil that acts as a reservoir for bacteria with the potential to thrive in freshwater environments (Crump et al., 2012). Assuming a low abundance of many of these taxa, both priority effects and ecological drift would arguably add to community stochasticity (Nemergut et al., 2013). By contrast, bacterial taxa occurring in 10–15 of the 20 ponds seemed to be assembled deterministically (Fig. 6A). Many dominant taxa of aquatic bacteria are relatively invariant at regional scales (Ostman et al., 2010). Dispersal by mass effect dynamics of abundant taxa (Urban, 2004) in combination with either an euryoecious or a stenoecious life style could thus lead to such deterministic patterns. The apparent return to stochasticity for OTUs occurring in almost all systems (Fig. 6A) may reflect a limitation of null model approaches if performed on a limited number of taxa (Chase et al., 2011). Only two bacterial OTUs were present in all ponds, one of them (Polynucleobacter, PnecC) being a well‐studied case of ecological diversification below the resolution of the rRNA gene (Hahn, 2003).

The changes in assembly patterns of the eukaryotic microbial communities at increasing frequency of occurrence did not show the same general trend as of bacteria (Fig. 6B), suggesting that fundamentally different mechanisms acted on rare and common bacteria and eukaryotic microbes. Eukaryotes in the surrounding soil might be less prone to colonize the ponds, as was reported for rivers and lakes (Crump et al., 2012). The most widely distributed OTU was affiliated with the cosmopolitan flagellate Cryptomonas marssonnii (Menezes and Novarino, 2003), more abundant in the turbid trough ponds (Mann–Whitney, W = 89, p = 0.005) and possibly favoured by its mixotrophic lifestyle (Saad et al., 2013).

GHG concentrations

Unlike many arctic ponds and lakes (Sobek et al., 2005; Repo et al., 2007; Breton et al., 2009), the majority of the studied ponds, in particular the polygonal ones, were net autotrophic, i.e. their water column CO2 levels were below the equilibrium concentration. While the quantity of Chl‐a was similar between pond types, the water transparency of the polygonal ponds (Fig. 1) allowed for formation of epibenthic cyanobacterial mats (A.L.M., personal observation). Such mats have been reported to be responsible for the bulk of primary production in comparable systems (Bonilla et al., 2005) and might explain the observed CO2 depletion. On the other hand, differences in Chl‐a were insufficient to explain why half of the trough ponds were CO2 undersaturated, whereas the other half was oversaturated. All these ponds were highly turbid, likely preventing cyanobacterial mats to form. Alternatively, differences in CO2 saturation levels could be related to chemical weathering of carbonates (Pierre et al., 2019). In contrast to previous reports (Breton et al., 2009), dissolved CO2 concentrations could not be predicted from the amount of DOC, suggesting that net auto‐ or heterotrophy was controlled by carbonate weathering or biotic factors such as food‐web structure (Schindler et al., 1997). The wide range of CO2 saturation levels might also be the result of asynchronous temporal dynamics, i.e. alternating blooms of autotrophs and grazers might cause ponds to switch between periods of net autotrophy and heterotrophy (Laas et al., 2012). While a conspicuous bloom of Chlamydomonas in a single trough pond (B05) paralleled by pronounced CO2 depletion seems to support this notion, a time‐resolved monitoring of CO2 concentrations and food web structure in several ponds would be required to substantiate it.

All ponds were oversaturated with CH4, highlighting the importance of thermokarst ponds as GHG sources (Laurion et al., 2010; Negandhi et al., 2014). The positive relationship between the isotopic signature of methane and the estimated abundances of methanotrophs reflects the accumulation of methane with higher 13C content in ponds with more methanotrophs. Such an accumulation is indicative of methanotrophy, as methanotrophs preferentially incorporate light methane (Smith et al., 1991). Methanotrophy is also suggested by the decrease of CH4 concentrations with the estimated abundance of methanotrophs (Fig. 4). Our findings expand previous results where a similar relationship was described at the sediment – water column interface (Crevecoeur et al., 2017) by highlighting the potential role of pelagic methanotrophs in mitigating CH4 emissions. However, methanotrophs featured a sizable community only in polygonal ponds, with an exceptionally high contribution to total bacterial reads in one pond (A07, 9%).

In contrast to prior findings from thermokarst ponds (Crevecoeur et al., 2017), the dominant methanotrophs in the polygonal ponds were of Type II. Type I and type II methanotrophs are distinct in their methanotrophic metabolic pathway as well as in their phylogenetic affiliation (Lüke and Frenzel, 2011). Type II methanotrophs have been reported to be dominant in boreal wetlands (Dedysh, 2009) and in some western Siberian lakes (Osudar et al., 2016). Although type II methanotrophs tend to thrive at higher temperatures than type I (Graef et al., 2011; He et al., 2012), the factors controlling the dominance of one type over the other are not fully understood (Osudar et al., 2016).

The low contribution of methanotrophs to total bacterial read numbers in the trough ponds illustrates that CH4 concentration from particular pond type may be largely unaffected by microbial consumption processes in the water column. CH4 fluxes in freshwater systems are controlled by both, biological and physical processes, such as temperature and stratification (Yvon‐Durocher et al., 2014; Bartosiewicz et al., 2015; Bartosiewicz et al., 2016). The relationship between the CH4 concentrations and δ2H‐H2O suggests that physical processes, such as evaporation losses could affect the potential release of CH4 in our study systems.

Conclusion

Shallow ponds across the Arctic Siberian tundra showed striking variability in their geochemical features and in their bacterial and eukaryotic microbial communities. While the geomorphological origin, habitat conditions and geographic location of these ponds clearly affected the occurrence of eukaryotic microbes as well as of the widely distributed bacterial taxa and methanotrophs, stochastic assembly processes shaped communities of ponds from the same type and geographic region. This stochasticity of assembly processes was reflected in the functional variability of the ponds, notably in their balance of net auto‐ versus heterotrophy and CH4 concentration. The overall relationship between planktonic methanotrophs and methane concentrations highlights the need for a better understanding of the habitat preference of methanotrophs to assess and predict carbon budgets in permafrost areas.

Experimental procedures

Study site

The field campaign took place from July 13 to 23 2018 at the Kytalyk Nature Reserve (70.82N, 147.47E), a low‐arctic tundra area in North‐East Siberia, Russian Federation. The research site within the reserve is located in the continuous permafrost region, in the Indigirka lowlands, 30 km north from Chokurdakh (Fig. 1). An extensive site description can be found in Juszak et al. 2016. The permafrost of the region is carbon‐ and ice‐rich, (Yedoma) (Grosse et al., 2013; Strauss et al., 2013). In total, 20 shallow ponds were sampled from two distinct geomorphological sites situated 1.5 km apart. At site A, we selected 10 ponds presumably formed by thermal ground contraction, creating low‐centred polygons that can be filled by precipitation (Lachenbruch, 1962) (subsequently referred to as polygonal ponds). At site B, we selected 10 ponds most likely formed by the thawing of the underlying permafrost ice, named trough ponds hereafter. Both ponds types are frequent in the region (Li et al., 2017). The ponds at each site, while located relatively close to each other, were selected according to their clear shore delimitation and their lack of direct hydrological connectivity. All ponds were sampled once over a period of 2 weeks. Pond B01 was sampled three times during the campaign to assess the variability of its physico‐chemical and biological properties (Fig. S6).

Physico‐chemical characterization of the ponds

The elevation of each pond was recorded by a hand‐held GPS device. The surface area was digitized with the ArcGIS software from an orthomosaic of drone images with a resolution of 15 cm taken during the sampling period. Maximum depth and active layer thickness in and around each pond were determined with a graduated steel rod. pH was assessed with pH paper (Whatman, indicator paper 0–14) and temperature with an analogue thermometer. Duplicate surface water samples were collected from all ponds with a Van Dorn sampler and filtered through pre‐combusted (4 h at 550°C) glass fibres filters (Whatman GF/F, pore‐size: 0.7 μm). The filtrate was collected in pre‐combusted glass vials and used to measure DOC, dissolved inorganic carbon (DIC), dissolved nitrogen (DN), optical properties of the dissolved organic matter (DOM) and the stable isotopic composition of the water (δ18O and δ2H). The filters served to determine particulate organic carbon (POC) and its stable isotopic composition (δ13C). Filters and filtrates were stored frozen in the permafrost cellar (≤ −5°C) of the research station and remained frozen until analysis in the laboratory. DOC, DIC and DN were measured using a Formacs HT‐I TOC/TN analyser (Skalar) calibrated with potassium phthalate. Optical properties of DOM were determined by spectrophotometric scans from 200 to 800 nm at steps of 1 nm (UV 1601, Shimadzu). The chromophoric DOM fraction (cDOM) was quantified by measuring the absorption coefficient of the water at 320 nm (a320). The specific ultraviolet absorbance index (SUVA254) was used as a proxy for the proportion of aromatic compounds (Weishaar et al., 2003). The Slope Ratio is an indicator of the molecular weight of the cDOM and was calculated across the spectral bands 275–295 and 350–400 nm as described in (Helms et al., 2008). To determine the stable isotope composition of the water, prefiltered samples were filtered again through 0.1 um pore size PES filters (Yeti, Infochroma) and measured by Off‐Axis Integrated Cavity Output Spectroscopy (OA‐ICOS) using a Liquid Water Isotope Analyser (LWIA‐24‐EP, Los Gatos Research). To determine POC and δ13C‐POC, the filters were dried for 72 h at 65°C, weighted and cut in half, then combusted at 900°C and analysed by Combustion Module‐Cavity Ring Down Spectroscopy (Picarro Inc.). The POC concentration was weighted by the filtered volume, ranging between 100 and 250 ml. Dissolved CO2 and CH4 concentration were obtained with the headspace technique (Kling et al., 1992; Laurion et al., 2010; Dinsmore et al., 2013). Fourty millilitre of pond water was equilibrated with 20 ml of ambient air. Samples were shaken, equilibrated for 1 min and the headspace was sampled in an evacuated helium‐purged 20 ml glass vial. The headspace and ambient air were collected in pentaplicates for each pond. CO2 and CH4 concentrations were measured by a gas chromatograph equipped with an FID detector and a methanizer and recalculated to aqueous values according to Henry's law:where KH is the Henry's constant adjusted for the water temperature and pGas the partial pressure of the gas in the headspace. Because samples were equilibrated in air, δ13C‐CH4 signatures were corrected for minor atmospheric contribution and further analyses were carried out only for samples with sufficiently high concentrations (CH4 > 0.2 μM), which resulted in missing information for two polygonal ponds. This analysis was performed in a pre‐concentration system (PreCon2018, Thermoscientific) coupled to an isotope ratio mass spectrometer (Gasbench II – Delta V Advantage, Thermoscientific). All environmental data were deposited at the Knowledge Network for Biocomplexity (doi:10.5063/F1P26WH3).

Biological components

The vegetation surrounding each ponds was assessed by determining the percent cover of each vascular plants and moss using a gridded quadrat of 50 × 50 cm2. Five quadrats were evaluated around each pond and the mean abundance of each species was used for subsequent analysis. Phytoplankton was collected in duplicates on glass fibre filters (Whatman GF/F, pore‐size: 0.7 μm), and stored frozen. Chl‐a was sequentially extracted with the following solvents: acetone: water (85:15 v/v), pure acetone and isopropanol: hexane (50:50 v/v). The extracts were vacuum concentrated, followed by volume reduction under a gentle stream of nitrogen, and stored at −80°C until their quantification with a U‐HPLC (Agilent 1290 Infinity) coupled with a diode array detector (Milani et al., 2019). Absorbance at 665 nm was used to measure Chl‐a concentration, identified with a standard (Sigma‐Aldrich). To estimate the abundance of the planktonic bacteria, unfiltered pond water was fixed with formaldehyde (2% v/v final concentration) and stored at 4°C for approximately 2–3 weeks before further analysis. Portions of 2 ml of the fixed water were filtered onto black polycarbonate filters (pore size: 0.22 μm) and stained with 4,6‐diamidino‐2‐phenylindole. A minimum of 800 cells per sample were counted under UV excitation by epifluorescence microscopy (Zeiss Axio Imager M1) at 1600× magnification (Posch et al., 2001).

In order to determine the composition of pond water microbial communities, 150 to 300 ml of pond surface water was filtered in duplicates onto 0.22 μm pore sized polyethersulfone filters (Millipore, GPWP). The filters were stored frozen until their analysis. The DNA was extracted with the DNeasy PowerBiofilm Kit (Qiagen, Germany) largely according to the manufacturer's specification, but with extending the inhibitor removal step to 1 h. The purified DNA was stored at −20°C in 10 mM Tris buffer. Partial ribosomal genes were amplified from the extracted DNA by PCR using the primer pair 799F‐1115R targeting bacteria with high specificity but excluding cyanobacteria and chloroplasts (Chelius and Triplett, 2001; Redford et al., 2010), and TAReuk454FWD1 – TAReukREV3 for eukaryotes (Stoeck et al., 2010). PCR products were sequenced on the Illumina MiSeq platform (2 x 300 bp, LGC genomics, Germany). Raw sequence data were deposited at the European Nucleotide Archive (PRJEB40506).

Raw data from rRNA gene sequencing were processed by an in‐house pipeline written in DELPHI, which included the joining of paired reads, quality control (0.5 average error rate per sequence), dereplication, exclusion of singleton reads, and a probabilistic clustering of operational taxonomic units (OTUs) comparable to UNOISE (Edgar and Flyvbjerg, 2015). Our algorithm first sorts sequences according to abundance, then assesses the probability that each queried sequence is derived from a more abundant ‘parent’, using published data to model the error probability of PCR amplification (propagating error) and sequencing (non‐propagating error). Then it assesses if the query sequence is chimeric using the whole data set (internal chimera check), and recursively assigns all sequences to their most likely ‘parents’. The remaining ‘ancestor’ sequences are retained as the most probable ‘true’ genotypes. Thus, there is no fixed distance threshold; instead, depending on the available information, sequence types with 2–3 mismatches can be discriminated. The pipeline was benchmarked using the Human Microbiome Mock Community (http://hmpdacc.org/HMMC/) against Uparse, Mothur, Qiime, Dada2 and Swarm, where it performed as well as Uparse in terms of accuracy.

The taxonomic assignment of OTUs was realized by a BLAST analysis of their representative sequences against the SILVA database (version 132) (Quast et al., 2012). Only bacterial OTUs with a query coverage ≥99% (i.e. 99% of bases from the OTU sequence could be aligned to a database sequence by BLAST) and eukaryotic OTUs with a query cover of 100% were included. Mitochondria and chloroplasts were removed from the bacterial OTU list, and spermatophytes and metazoans from the eukaryotic OTU list. A total of 3156 and 2968 OTUs were detected prior to normalization for bacteria and eukaryotes respectively. The OTU lists were deposited at the Knowledge Network for Biocomplexity (doi:10.5063/F1P26WH3).The two datasets were normalized to equal sample size using the ‘rarefy_even_depth’ function (R package ‘Phyloseq’) that randomly removes reads from the dataset using replacement (McMurdie and Holmes, 2013). This resulted in 2466 bacterial and 1591 eukaryotic OTUs with a sequencing depth of 7822 and 1460 reads per sample, respectively. The estimated abundances of methanotrophs were calculated from the proportions of their read numbers multiplied by total bacterial cell numbers as obtained from microscopic counts.

Data analysis

All statistical analyses were carried out with the R software (version 3.6.1). Differences in environmental parameters between the two sets of ponds were tested by Mann–Whitney or t‐tests, depending on the normality and homoscedasticity of the data. False discovery rate of multiple testing was corrected using the Benjamini‐Hochberg method (Benjamini and Hochberg, 1995).

The β‐diversity of the communities was assessed using the BC similarity index, defined as 1 – BC dissimilarity with the ‘Vegan’ package (Oksanen et al., 2013). Duplicate microbial communities from single ponds were analysed separately, resulting in 20 data points per sites. BC similarity of bacterial and eukaryotic communities in such duplicate samples was high (mean BC similarity in polygonal ponds: bacteria 81%; eukaryotes 81%; in trough ponds: bacteria 71%; eukaryotes 76%). Moreover, samples taken from the same ponds were significantly more similar to each other than to samples taken from different ponds (Chi‐square tests: bacteria: χ 2 = 82, p < 0.01; eukaryotes: χ 2 = 115, p < 0.01).

Non‐metric multidimensional scaling was performed using ‘metaMDS’ and PERMANOVAs were computed with the ‘adonis’ function from the ‘Vegan’ package. The relative contribution of OTUs from different phylogenetic lineages to the average dissimilarity of polygonal and trough pond communities was assessed with the ‘simper’ function from the same package (Clarke, 1993). One replicate sample of pond B09 only contained a single eukaryotic read and was, therefore, excluded from the analysis. Similarity profile analysis (‘simprof’ with 999 permutations) from the ‘Clustsig’ package was used to assess the significance of dissimilarity between the samples (Clarke et al., 2008). The relationship between the collected environmental variables and the microbial communities was inferred by Mantel tests on the duplicates‐pooled datasets. A db‐RDA was computed on the pooled dataset with the function ‘capscale’ from the ‘vegan’ package to analyse which physico‐chemical parameters were the most important in shaping microbial community structure. The model was built on a reduced number of variables chosen with the function ‘bioenv’. Pond B05 was excluded due to its unique phytoplankton bloom.

In order to investigate the assembly processes of the microbial communities, the modified RC index that is based on read numbers (Chase et al., 2011) and the NST developed by Ning et al., (2019) were calculated using the ‘NST’ package. These indices use null model‐based approaches to quantify ecological stochasticity. The RC metric takes into account variability of α‐diversity; it indicates if two communities are more (−1 < RC < −0.95) or less (0.95 < RC < 1) similar than expected by chance or if stochastic processes dominate (−0.95 < RC < 0.95). Similarly, if NST < 50%, community assembly is predominantly governed by deterministic processes, while NST values >50% indicate the dominance of stochastic processes. For these analyses, the replicates from each pond were pooled. The NST index was also calculated separately for subcommunities of OTUS that were present in an increasing number of systems. This allows to assess the relative importance of stochastic and deterministic processes with respect to rarely occurring and more ubiquitous community members.

We used the βNTI to assess phylogenetic turnover across the ponds (Stegen et al., 2012; Stegen et al., 2013). For this purpose, the bacterial and eukaryotic OTUs were first independently aligned with PASTA using default settings (Mirarab et al., 2014), and maximum‐likelihood phylogenies were constructed using IQ‐TREE (Nguyen et al., 2015) with the best‐fitting substitution models chosen by ModelFinder (Kalyaanamoorthy et al., 2017) based on the Bayesian Information Criterion. The βNTIs were then calculated from phylogenetic distances using the code provided in (Stegen et al., 2013) with the ‘picante’ package. Communities are more closely related than expected by chance when βNTI >2, less closely related when βNTI < −2 and randomly related when −2 < βNTI < 2. To understand the impact of evaporation on the methane dynamics, correlation between δ2H‐H2O and methane concentration as well as δ13C‐CH4 was performed. The role of methanotrophs was assessed with correlations between the estimated abundance of methanotrophs and methane concentration as well as δ13C‐CH4. Additionally, methane concentrations and estimated abundance of methanotrophs were Box‐Cox transformed (Box and Cox, 1964) and a regression of methane concentration on estimated abundance of methanotrophs was calculated. Pond B05 was excluded from this analysis because there were no reads from methanotrophs in this system.

Fig. S1. Orthomosaic of drone images of (A) site A, (B) site B. Points and label represent individual ponds sampled in this study.

Fig. S2. Hierarchical agglomerative clustering of the ponds based on the surrounding vegetation (distance: Bray‐Curtis, agglomeration: average). Colours represent significant clusters detected by similarity profile analysis (Simprof).

Fig. S3. Isotopic signature of pond water shown in relation to the Global Meteoric Water Line (GMWL).

Fig. S4. Phylogenetic relationship of methanotrophs (Maximum Likelihood method) and their relative abundance in each pond. Each circle represents one pond.

Fig. S5. db‐RDA of the (A) bacterial (B) micro‐eukaryotic communities. Shaded points represent the OTUs and arrows represent significant environmental parameters. On the axis, the percentage of inertia is represented for the fitted model and for the total inertia. Pond B05 was excluded from the analysis of the bacterial communities.

Fig. S6. A. Coefficient of variation of the physico‐chemical properties in pond B01 from the three sampling dates: t1: 13/07/2018, t2: 17/07/2018 and t3: 23/07/2018. Abbreviations are described in Table 1.

B. Bray‐Curtis dissimilarity for the prokaryotic and eukaryotic communities for the three sampling dates.

Table S1. Average read numbers of OTUs that significantly contributed to the difference of the bacterial microbial communities between polygonal and trough ponds (SIMPER analysis). Only OTUs with a contribution >1% are shown. Distance: percentage difference to closest know sequence (accession number).

Table S2. Average read numbers of OTUs that significantly contributed to the difference of the eukaryotic microbial communities between polygonal and trough ponds (SIMPER analysis). Only OTUs with a contribution >1% are shown. Distance: percentage difference to closest know sequence (accession number).

Click here for additional data file.