Expansion of shrubs could result in local loss of soil bacterial richness in Western Greenland

ABSTRACT

Climate warming in Greenland is facilitating the expansion of shrubs across wide areas of tundra. Given the close association between plants and soil microorganisms and the important role of soil bacteria in ecosystem functioning, it is of utmost importance to characterize microbial communities of arctic soil habitats and assess the influence of plant edaphic factors on their composition. We used 16S rRNA gene amplicons to explore the bacterial assemblages of three different soil habitats representative of a plant coverage gradient: bare ground, biological soil crusts dominated by mosses and lichens and vascular vegetation dominated by shrubs. We investigated how bacterial richness and community composition were affected by the vegetation coverage, and soil pH, moisture and carbon (C), nitrogen (N) and phosphorus (P) contents. Bacterial richness did not correlate with plant coverage complexity, while community structure varied between habitats. Edaphic variables affected both the taxonomic richness and community composition. The high number of Amplicon Sequence Variants (ASVs) indicators of bare ground plots suggests a risk of local bacterial diversity loss due to expansion of vascular vegetation.

Arctic, environmental filtering, climate change, edaphic parameters, metabarcoding, 16S rRNA gene

INTRODUCTION

In many Arctic environments, rising temperatures are causing an increased melting of glacial ice and an expansion of shrubs (Sturm, Racine and Tape 2001; Tape, Sturm and Racine 2006) at the expense of vegetation dominated by bryophytes and lichens (Normand et al. 2013; Cahoon, Sullivan and Post 2016; Vowles and Björk 2019). These trends are particularly apparent in Greenland, one of the regions in the Arctic most affected by climate change (Howat and Eddy 2011; Bevis et al. 2019).

Soil microbial communities and vegetation are linked by a variety of direct and indirect interaction: plants provide photosynthetically fixed carbon and low molecular weight root exudates that are used as energy sources by soil microorganisms, whose community composition is in turn shaped by their ability to metabolize different compounds and to resist different antimicrobial metabolites present in the exudates (Marschner, Crowley and Yang 2004; Berg and Smalla 2009; Berendsen, Pieterse and Bakker 2012; Schulz-Bohm et al. 2018). On the other hand, rhizosphere microorganisms can strongly influence plant growth and health. Arctic tundra is one of the most nitrogen limited environments on Earth and the input from biological nitrogen fixation, mainly contributed by Cyanobacteria, is critical to the development of these ecosystems (Solheim et al. 2006). Additionally, the mobilization of nitrogen from organic matter is an important source for plant growth, and is at least partly influenced by bacterial activity (Chen et al. 2014; Leff et al. 2015). Furthermore, plant growth promoting bacteria may improve resource acquisition, modulate plant hormones secretion or act antagonistically against pathogens (Glick 2012; Berendsen, Pieterse and Bakker 2012). This close association between plants and microorganisms have led to the concept of plants as ‘superorganisms’ that partly rely on their microbiome interactions for specific functions and traits (Mendes, Garbeva and Raaijmakers 2013).

The reported decline of mosses and lichens due to climate change and their replacement by vascular plants, particularly shrubs (Normand et al. 2013; Cahoon, Sullivan and Post 2016; Vowles and Björk 2019), can be expected to affect soil bacteria. Changes in the composition and structure of soil bacterial communities in the Arctic have been recently documented in simulations of long-term warming (Deslippe et al. 2012), higher nutrient availability (Koyama et al. 2014; Männistö et al. 2016) and altered precipitation regimes (Ricketts et al. 2016), as well as in natural successions in glacier forelands (Kwon et al. 2015; Kim et al. 2017) and along permafrost thaw gradients (Deng et al. 2015; Frank-Fahle et al. 2014). These responses of bacterial communities to climate changes, and their impact on ecosystem functionality are of utmost scientific relevance. Despite this importance, our understanding of the warming-driven changes in Arctic soil habitats remains uncomplete. For example, the potential effects of the shrub expansion on the structure and functioning of associated bacterial communities are challenging to accurately forecast, since the effects of vegetation type on soil communities in Arctic tundra are not well known (Krab et al. 2019). Differences have been observed across Alaskan tundra vegetation types (Wallenstein, McMahon and Schimel 2007), suggesting that plant communities influence bacterial communities via the quantity and quality of the litter supply, and by modifying the soil physical environment. In soils of Eastern Greenland (Ganzert, Bajerski and Wagner 2014) abiotic parameters, related to different habitats, shaped microbial communities. In a Canadian low Arctic tundra system Denaturing Gradient Gel Electrophoresis analyses revealed that vegetation coverage plays a key role in shifting bacterial communities (Chu et al. 2011). A more recent study in the Canadian Arctic, using 16S rRNA gene pyrosequencing, suggested that soil responses to warming would be vegetation-specific, likely due to the differences in the structure of microbial communities associated with different plants (Shi et al. 2015). However, other works in Finnish (Männistö, Tiirola and Häggblom 2007) and Canadian tundra ecosystems (Buckeridge et al. 2010), did not show any differences in bacterial communities among vegetation types.

Given the wide metabolic capabilities of bacteria, a deeper knowledge of the relationships between biotic and abiotic factors shaping soil bacterial community structure and function may be useful to predict the effects of global change on the vast and highly vulnerable Arctic soil carbon stocks (Crowther et al. 2016). Shrub encroachment, especially of deciduous species, such as Betula nana, may produce more labile compounds from leaf litter, increasing the turnover of soil C (Weintraub and Schimel 2005; Wookey et al. 2009). The impact of warming will depend on how efficiently plant-derived carbon is incorporated into microbial biomass or converted to carbon dioxide and released to the atmosphere (Cotrufo et al. 2013). In addition, plant litter and root exudates can facilitate the activity of microorganisms with enhanced decomposition abilities for old stocks of organic matter, in a process known as priming (Fontaine, Mariotti and Abbadie 2003; Kuzyakov 2002; Walker et al. 2015).

In our study, we selected three different habitats in Western Greenland, representing a gradient of vegetation complexity: bare ground (BG), biological soil crusts (BSCs) dominated by mosses and lichens and vascular vegetation (VV) dominated by shrubs, e.g. Empetrum, Vaccinium, Betula and Salix. The main aims were: (i) to characterize the diversity and composition of soil bacterial communities in this Arctic region; (ii) determine whether the diversity and community composition of soil bacteria are related to edaphic parameters and (iii) identify the abiotic parameters most closely associated with variation in bacterial communities. Additionally, focusing only on the vascular vegetation plots, we aimed to test for correlations between relative cover of different shrub genera and the composition of soil bacterial communities, to assess the impact of shrub expansion on these patterns.

MATERIALS AND METHODS

Sampling

Sampling was carried out in July 27–31, 2017 in the area of Kobbefjord, Nuuk, West Greenland (64°08’ N, 51°23’ W). The climate of the area is classified as low Arctic (Jonasson et al. 2000). The mean annual air temperature in the years 2008–2010 was 0.7°C, with the mean air temperature of the warmest month, July, 10.7°C. Over the same period annual precipitation ranged between 838 and 1127 mm, with an average of 25–50% of the total annual precipitation falling as snow during the winter period (Søndergaard et al. 2012). Samples were collected in the area close to the NERO line, a permanent vegetation transect established in 2007 to monitor changes in vegetation species composition (Bay, Aastrup and Nymand 2008).

In total, 20 plots (2m²), representing three habitat types, were sampled: five in bare ground (BG), six in biological soil crusts (BSC) and nine in soils covered with vascular vegetation (VV). Exact coordinates, elevation and shrub genera composition for each plot are listed in Table 1. In each plot, three replicate soil samples (up to 10 cm depth) were collected aseptically, after removing the top of the soil, plant litter in vascular vegetated plots and the superficial coverage of mosses and lichens in BSC plots. Samples were transported in sterile bags and stored at −20°C at University of Tuscia, Italy, until further processing.

Edaphic parameters

Gravimetric soil water content was measured on 5 g subsamples dried at 105°C (Reynolds 1970). Measurements were repeated until no variation in weight was observed. pH was measured in a 1:2.5 suspension of dried soil in deionized water, with a HI9321 pH meter (Hanna Instruments Woonsocket, RI). For each sample, soil moisture and pH were measured in independent triplicates. Phosphorus (P), carbon (C) and nitrogen (N) content were analysed at Eger Innovations Nonprofit Kft. (Eszterházy Károly University, Eger, Hungary). P content was measured by Microwave Plasma Atomic Emission Spectrometry (MP-AES) and C and N content by CNS elemental analyser.

DNA extraction, amplification and sequencing

For each sample, DNA was extracted from 0.5 g of soil using DNEasy Powersoil kit (QIAGEN, Hilden, Germany), according to the manufacturer's protocol. The V4 hypervariable region of the 16S rRNA gene was amplified using 515F (Parada, Needham and Fuhrman 2016) and 806R (Apprill et al. 2015) primers; libraries were prepared following the protocol of Minich et al. (2018). The equimolar pool of uniquely barcoded amplicons was paired-end sequenced (2 × 300 bp) on an Illumina MiSeq platform at the Vincent J. Coates Genomics Sequencing Laboratory at University of California, Berkeley.

Bioinformatic analyses

Bcl files were converted to Fastq files, demultiplexed and primer removed using bcl2fastq (v 2.18). Dual-matched 8-bp indexes were used to eliminate the occurrence of “barcode bleed” (or tag-switching) between samples.

Demultiplexed 16S rRNA gene sequences were processed with the QIIME2 (Quantitative Insights Into Microbial Ecology, v. 2018.11; Bolyen et al. 2018) platform. 3179 253 starting sequences were denoised, trimmed to length 160 bp, merged and clustered in Amplicon Sequence Variants (ASVs), using DADA2 (Callahan et al. 2016), which includes phiX reads removal and chimera detection. We obtained 12 143 quality filtered ASVs, each with at least two reads in the total dataset. Taxonomy was assigned with the q2-feature-classifier within the database Greengenes v. 13_8 (99% OTUs from 515F/806R region of sequences). Chloroplasts, mitochondrial, chimeric and low identity ASVs (less than 80% identity to other prokaryotic 16S rRNA sequences) were removed, retaining 10 980 ASVs. The dataset was normalized for subsequent analyses, rarefying the number of reads per sample to the lowest reads obtained (14 251 reads) using the rrarefy function in the vegan package v. 2.5-2 (Oksanen et al. 2018) in R v. 3.5.2 (R Core Team, 2018), retaining a total of 10 578 ASVs. Sequences of ASVs were submitted to NCBI gene bank (BioProject PRJNA550020).

Statistical analyses

Unless otherwise specified, all statistical analyses were carried out with the vegan package v. 2.5-2 (Oksanen et al. 2018) in R v. 3.5.2 (R Core Team, 2018). Total bacterial richness (including all the ASVs retrieved), as well as relative richness (proportion of ASVs in a sample belonging to each group) and relative abundance (proportion of total reads in each sample assigned to each group) of most abundant phyla and classes among the three habitats were compared using ANOVA and Tukey's HSD test. Linear regression analyses were used to examine relationships between edaphic factors (pH, soil moisture, C, N and P content, and C/N ratio) and bacterial phyla and classes relative richness and abundances

We performed Non-Metric Multidimensional Scaling (NMDS) of the weighted Bray-Curtis distances of Hellinger transformed matrix of the bacterial community. We used the envfit R function to project edaphic variables (pH, soil moisture, C, N and P content, and C/N ratio) and the relative abundance values of the shrub genera (Betula, Empetrum, Salix and Vaccinium) or of different taxonomic groups onto the NMDS ordinations. In addition, we tested whether bacterial communities were statistically different among habitat types using the multi response permutation procedure (MRPP). We determined preferences of unique ASVs for each habitat using indicator species analyses on the Hellinger transformed matrix of the bacterial community (Dufrêne and Legendre 1997) in PC-ORD v. 6.0 (McCune et al. 2002).

Permutational multivariate analysis of variance (PerMANOVA; Anderson 2001) was carried out on Hellinger transformed Bray-Curtis distance matrix to determine the effect of each soil physicochemical parameter on the observed variance of the total community and of dominant phyla. Significant variables obtained from this analysis, were considered in a model to determine the combined effect of soil parameters on the variance of the community. The same approach was used taking into account only VV plots (27 samples in total) in order to assess the effect of the relative abundance of the four dominant shrub genera (Salix, Betula, Vaccinium and Empetrum), in combination with edaphic parameters, on the variance of the total community and of the different phyla.

RESULTS

Bacterial richness and abundance patterns

The filtered and rarefied dataset of the 16S rRNA gene amplicons contained 10 578 bacterial ASVs. The proportions of ASVs found exclusively in VV plots was the highest (41.6%) compared to BSC and BG plots (12.1 and 18.9%, respectively), whereas the BG samples showed the highest number of indicator ASVs (407 ASVs, compared with 86 and 114 for VV and BSC samples, respectively). Only 759 ASVs (7.2%) out of 10 578 were present in all the habitats.

The total richness of the bacterial communities was higher in VV and BG samples compared to BSC, with the two former not significantly different from each other (Fig. 1a). We selected 10 most abundant phyla of 35 identified (encompassing both Bacteria and Archaea) and 12 most abundant classes of 103 retrieved for further analyses. The 10 phyla selected represented more than 80% of total reads identified and more than 95% in many samples. The 12 classes were representative of more than 65% of reads identified at this level in all the samples.

Figure 1.

Richness of the total bacterial community and relative richness of the 12 dominant bacterial classes in each habitat (green, Bare Grounds plots; red, Biological Soil Crusts plots; blue, Vascular Vegetation plots). Letters indicate significant differences in one-way ANOVA post-hoc Tukey's HSD test (P <0.05).

Among the most representative classes and phyla studied, those that had the highest richness in BG plots, were Ktedonobacteria (Chloroflexi), Spartobacteria (no statistical significance with BSC; Verrucomicrobia) and the phylum Cyanobacteria (Fig. 1; Figure S1, Supporting Information).

Classes with the highest richness in BSC plots were Acidobacteriia and Solibacterales (Acidobacteria), while DA052 within the same phylum showed the lowest richness and abundance in this habitat. Classes Actinobacteria and Thermoleophilia (Actinobacteria) had higher richness and abundance in BSC plots, compared to BG and VV. The same was for the class Alphaproteobacteria and the phylum Armatimonadetes, that had higher richness for this habitat type. Instead, Planctomycetia (Planctomycetes), the class Betaproteobacteria and the phylum Chlamydiae showed lower richness in BSC plots (Fig. 1; Figure S1, Supporting Information).

Gammaproteobacteria and Deltaproteobacteria (Proteobacteria), as well as Bacteroidetes had the highest richness in VV plots (Fig. 1; Figure S1, Supporting Information). The same general trends were observed in the relative abundances of individual taxonomical groups (Figures S2 and S3, Supporting Information).

Effect of environmental parameters on richness and abundance of different groups

Soil moisture, P and N content of soil samples showed an increase from BG plots to vegetated plots (BSC + VV), while pH decreased. C content and the C/N ratio also increased from BG plots to the vegetated plots, but with higher values in BSC compared to VV plots (Fig. 2). C content and the C/N ratio were the main predictors of total bacterial richness, with a negative correlation (slope = −4.06 and r² = 0.105 for C content; slope = −10,01 and r² = 0.160 for C/N ratio), while N content was only marginally significant (slope = −72.89 and r² = 0.035), and both soil moisture and pH were positively correlated (slope = 2.14 and r² = 0.052 for soil moisture; slope = 101.38 and r² = 0.058 for pH; Figure S4, Supporting Information).

Figure 2.

Differences in soil parameters across the three habitats (green, Bare Ground plots; red, Biological Soil Crusts plots; blue, Vascular Vegetation plots). Letters indicate significant differences in one-way ANOVA post-hoc Tukey's HSD test (P <0.05).

Richness and abundance of the dominant phyla and classes were significantly related to both C and N content and C/N ratio. Correlations were positive for both richness and abundance of Acidobacteria, Actinobacteria (N content not significant for richness), Armatimonadetes and Bacteroidetes and with the richness of Proteobacteria and Verrucomicrobia. Conversely, these three parameters were negatively correlated with both the richness and abundance of Chloroflexi, with the richness of Cyanobacteria and with the abundance of Planctomycetes and Verrucomicrobia (Fig. 3, Figure S5, Supporting Information). pH significantly influenced phyla richnesses and abundances. There were negative correlations with both the richness and the abundance of Acidobacteria, Actinobacteria, Armatimonadetes and Bacteroidetes and with the richness of Proteobacteria and Verrucomicrobia; and positive correlations with both richness and abundance of Chloroflexi, with the richness of Cyanobacteria and with the abundance of Verrucomicrobia (Fig. 3, Figure S5, Supporting Information). Soil moisture and P content had a significant effect on a smaller number of the dominant taxonomic groups. Soil moisture was positively correlated with the richness and abundance of Bacteroidetes, Chlamydiae and Proteobacteria and negatively with Chloroflexi (Fig. 3, Figure S5, Supporting Information). P content was positively correlated with both richness and abundance of Bacteroidetes, the richness of Proteobacteria and Verruomicrobia and the abundance of Actinobacteria. It was also negatively correlated with both richness and abundance of Chloroflexi, with the richness of Cyanobacteria and with the abundance of Verrucomicrobia (Fig. 3, Figure S5, Supporting Information). The trends were similar for the classes analysed within these phyla (see Table S2, Supporting Information for all statistical details).

Figure 3.

Summaries of linear regression models for the variation of richness of 9 dominant bacterial phyla in relation to soil parameters: soil moisture, pH, P content, C and N content and C/N ratio. The significance of the regressions is indicated as *** P < 0.001, ** P < 0.01, *P <0.05. Only significant regressions are reported. See Figure S4 (Supporting Information) for individual plots with data points shown.

Community composition

Bacterial communities structure was well differentiated between the habitat types (Fig. 4a, b; MRPP P = 0.001, A = 0.132). The same differentiation was also apparent when analysing the phyla-level composition (Figure S6, Supporting Information). Among variables fitted to the ordinations, all soil parameters, relative abundances of the four main shrub genera (Table S1, Supporting Information) and of the different phyla, were significant, except for the relative abundances of Betula, Empetrum and Cyanobacteria (Table S3, Supporting Information).

Figure 4.

Nonmetric multidimensional scaling (NMDS) ordinations of the differences (Bray–Curtis distance) in composition of bacterial communities (Hellinger transformed ASVs abundances) in the habitats studied (green, Bare Ground plots; red, Biological Soil Crusts plots; blue, Vascular Vegetation plots) for the total bacteria communities. Arrows represent projections of (a) edaphic variables (pH, soil moisture, carbon, nitrogen and phosphorus content, and C/N ratio) and the relative abundance values of the shrub genera (Betula, Empetrum, Salix and Vaccinium); and (b) relative abundances of dominant phyla.

When the effects of the habitat type (BG, BSC or VV) and single edaphic parameters were tested on both the total community composition and the phyla-level composition, all the variables were significant, with the habitat type always explaining the highest proportion of variation (26% for the total community; Table S4, Supporting Information). Conversely, when the variables were combined additionally in a model, only the type of habitat, the pH and the C/N ratio had an independent effect on the observed variance for both the total community and all the phyla considered, with the only exception of Armatimonadetes (Table 1). Both for the total communities and many phyla, N content resulted to be an independent parameter, whereas for other phyla (Bacteroidetes, Chlamydiae, Chloroflexi and Cyanobacteria) soil moisture was determinant. C and P content were never explaining an additional variance than the other parameters when combined additionally to them (Table 1).

Table 1.

Proportion of variation in bacterial community composition explained by habitat (categorical) and soil variables (continuous), based on permutational multivariate analyses of variance. Variables with significant results in individual analyses (Table S4, Supporting Information) were added sequentially in reverse order on explained variance to a combined model. Variables that remained significant in the combined model are in bold.

All Bacteria			Acidobacteria			Actinobacteria			Armatimonadetes
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Habitat	26.730	0.0001	Habitat	33.178	0.0001	Habitat	25.249	0.0001	Habitat	14.178	0.0001
pH	4.400	0.0001	pH	4.937	0.0001	pH	4.275	0.0009	pH	3.288	0.0004
C/N ratio	3.217	0.0027	C/N ratio	3.323	0.0030	C/N ratio	4.223	0.0006	C/N ratio	1.718	0.2851
C	1.439	0.2206	C	1.251	0.2634	C	1.799	0.1010	Soil moisture	1.643	0.3552
N	2.965	0.0068	N	2.789	0.0109	N	3.579	0.0019	C	1.311	0.7322
P	1.335	0.2886	P	1.197	0.2934	Soil moisture	1.594	0.1640	N	1.997	0.1165
Soil moisture	1.669	0.1225	Soil moisture	1.579	0.1255	P	1.340	0.2904	P	1.405	0.6362
Residuals	58.244		Residuals	51.747		Residuals	57.941		Residuals	74.460
Bacteroidetes			Chlamydiae			Chloroflexi			Cyanobacteria
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Habitat	21.732	0.0001	Habitat	6.754	0.0001	Habitat	20.075	0.0001	Habitat	9.266	0.0001
pH	4.253	0.0001	Soil moisture	2.946	0.0006	pH	3.496	0.0016	C/N ratio	2.819	0.0114
C/N ratio	2.456	0.0199	N	2.010	0.1504	C/N ratio	3.428	0.0015	pH	2.498	0.0370
Soil moisture	3.158	0.0032	pH	2.444	0.0100	C	1.536	0.2857	C	1.549	0.5352
C	1.537	0.2281	P	1.603	0.7579	Soil moisture	3.149	0.0034	Soil moisture	3.127	0.0027
N	1.849	0.1051	C	2.037	0.1321	N	1.790	0.1519	N	4.227	0.0002
P	1.487	0.2678	C/N ratio	2.389	0.0152	P	1.327	0.4641	P	1.510	0.5799
Residuals	63.527		Residuals	79.816		Residuals	65.199		Residuals	75.003
Planctomycetes			Proteobacteria			Verrucomicrobia
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Habitat	20.573	0.0001	Habitat	25.647	0.0001	Habitat	29.161	0.0001
C/N ratio	3.591	0.0010	pH	4.157	0.0003	pH	5.334	0.0002
pH	3.342	0.0027	C/N ratio	3.320	0.0027	C/N ratio	2.973	0.0003
C	1.499	0.2898	C	1.285	0.3321	C	1.734	0.0774
N	2.307	0.0338	N	3.379	0.0013	N	2.761	0.0103
P	1.163	0.6371	Soil moisture	1.800	0.0912	P	1.531	0.1374
Soil moisture	1.339	0.4134	P	1.431	0.2318	Soil moisture	1.506	0.1771
Residuals	66.187		Residuals	58.983		Residuals	55.001

Effect of shrub coverage on bacterial community composition

PerMANOVA analysis was used to test for the effect of shrub community composition on soil communities. When the abundance of different shrubs and the soil parameters were tested independently, Salix coverage explained the highest variance in community composition for many groups (Table S5, Supporting Information), and Betula coverage, never significant, was the main determinant only for Cyanobacteria community composition (Table 2). pH was the second strongest predictor, explaining the highest variance in community composition for many groups (Table 2). Finally, C, N and soil moisture, and Empetrum abundance had a significant independent effect on the variance of total community and, to a different extent, on the phyla studied (Table 2).

Table 2.

Proportion of variation in bacterial community composition in the vascular vegetation plots explained by soil variables and relative abundance of shrub genera, based on permutational multivariate analyses of variance. Variables with significant results in individual analyses (Table S6, Supporting Information) were added sequentially in reverse order on explained variance to a combined model. Variables that remained significant in the combined model are in bold.

All Bacteria			Acidobacteria			Actinobacteria			Armatimonadetes
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Salix	16.467	0.0001	Salix	21.028	0.0001	Salix	20.615	0.0001	Salix	9.416	0.0001
pH	7.193	0.0013	pH	7.184	0.0039	pH	7.667	0.0005	pH	4.565	0.2019
C	6.068	0.0073	C	6.422	0.0103	Vaccinium	9.157	0.0006	Residuals	86.019
Empetrum	4.867	0.0307	N	3.369	0.2002	Empetrum	4.512	0.0399
C/N ratio	3.305	0.2331	C/N ratio	3.901	0.1180	C	3.258	0.1665
N	4.606	0.0442	P	1.958	0.7156	N	1.904	0.6643
Vaccinium	2.490	0.5848	Vaccinium	3.192	0.2379	Soil moisture	7.241	0.0009
Soil moisture	4.644	0.0426	Empetrum	3.366	0.1984	C/N ratio	2.202	0.5044
P	2.839	0.4011	Soil moisture	4.891	0.0428	P	2.332	0.4502
Residuals	47.521		Residuals	44.689		Betula	2.561	0.3713
						Residuals	38.551
Bacteroidetes			Chlamydiae			Chloroflexi			Cyanobacteria
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Salix	17.137	0.0001	Empetrum	7.805	0.0001	Salix	12.335	0.0001	Betula	7.118	0.0091
pH	5.227	0.0566	Vaccinium	5.855	0.0034	Vaccinium	9.338	0.0003	Salix	6.257	0.154
Empetrum	4.988	0.0713	Salix	5.608	0.0045	Empetrum	5.505	0.0167	Residuals	86.626
C	4.594	0.1047	Soil moisture	4.516	0.0635	pH	4.155	0.1127
Residuals	68.053		N	3.531	0.4173	C	4.113	0.1166
			Residuals	72.684		Soil moisture	6.741	0.0020
						N	1.893	0.9105
						P	3.270	0.3373
						C/N ratio	1.946	0.9000
						Residuals	50.704
Planctomycetes			Proteobacteria			Verrucomicrobia
Variable	Variance (%)	p	Variable	Variance (%)	p	Variable	Variance (%)	p
Salix	12.656	0.0001	Salix	15.488	0.0001	pH	15.560	0.0001
pH	6.477	0.0046	pH	7.236	0.0007	Salix	8.908	0.0002
Vaccinium	7.521	0.0019	C	5.573	0.0113	C	7.659	0.0005
C	2.953	0.5087	Empetrum	5.308	0.0148	N	2.675	0.4706
C/N ratio	3.843	0.1747	C/N ratio	3.681	0.1367	Soil moisture	6.714	0.0029
Empetrum	3.938	0.1549	Vaccinium	3.729	0.1319	P	3.104	0.2894
N	3.536	0.2572	Soil moisture	5.806	0.0082	C/N ratio	2.549	0.5301
P	3.024	0.4711	N	2.584	0.5444	Empetrum	4.447	0.0503
Residuals	56.051		P	2.799	0.4284	Vaccinium	1.774	0.8979
			Residuals	47.797		Residuals	46.610

DISCUSSION

Our study reports the differences in richness and community composition along a gradient of vegetational complexity in Western Greenland. Our findings show that the structure of soil bacterial communities is strongly influenced by vegetation complexity and by the identity of dominant shrub genera. Our results partly confirm previously observed patterns, but also offers new insights. For example, Wallenstein, McMahon and Schimel (2007) found that bacterial communities in Alaskan acidic tundra soils were dominated by Acidobacteria in tussock tundra, while Proteobacteria dominated shrub tundra. In our study, Acidobacteria clearly preferred BSC soils that had the lowest pH, but they were much less diverse in vascular vegetation plots with similar pH that were dominated by dwarf shrubs (Figure S1, Supporting Information). For Proteobacteria, our data showed that phylum-level patterns can mask potentially important class-level differences. Out of the four Proteobacteria classes, three showed highest richness in shrub-dominated plots, while Alphaproteobacteria showed highest richness in BSC soils (Fig. 1), a pattern identical to that of Acidobacteria. Acidobacteria are generally considered k-strategists, with lower growth rates, but high efficiency in converting nutrients to biomass and high tolerance to toxic compounds. This results in a greater ability to compete in oligotrophic environments, which accords with their preference for BSC plots (Kielak et al. 2016). Instead, Proteobacteria are generally considered copiotrophic organisms. In the Arctic tundra, organisms of this phylum have been reported to be more abundant after fertilization experiments (Koyama et al. 2014), likely due to increased organic matter input by vascular plants (Ramirez et al. 2010).

The highest number of indicator ASVs recorded in BG plots suggests that these habitats harbour a unique pool of bacteria, adapted to thrive in these conditions. Similar patterns have been reported for fungi in Eastern Greenland (Grau et al. 2017). With the metabarcoding approach it is not possible to estimate the proportion of sequences derived from relic DNA, often considered to be the majority of the total DNA in low biomass soils (Carini et al. 2016). Although we cannot rule out the possibility of some of the observed patterns to be due to such relict sequences, the sequences found exclusively in BG plots in the dataset often belonged to stress-tolerant taxa, such as those of the phylum Chloroflexi. If these taxa are typical of bare-ground habitats, the recent expansion of shrubs into these habitats would likely result in the local extinction of these stress-tolerant microbes, since they could be expected to be outcompeted in fully vegetated habitats.

The highest richness observed for BG and VV plots relative to BSC (Fig. 1) differs from results found on a primary successional gradient of an Arctic glacier foreland, where richness positively correlated with vegetation complexity (Kwon et al. 2015). Our results also differs from previous comparisons of vegetated and non-vegetated soils, where lower diversity was observed in the former (Tam et al. 2001; Kumar et al. 2016). In this study, plant coverage is among the most crucial environmental factors influencing bacterial community composition.

In our study, a significant proportion of bacterial community variance was also explained by edaphic parameters (e.g. pH, C/N ratio), that were different among habitat types. Of these, soil pH was the best predictor of community composition and had a key role in predicting the richness and relative abundance of many taxonomic groups (Table 1). In agreement with other studies on arctic communities (Chu et al. 2010; 2011; Männistö, Tiirola and Häggblom 2007), the effect of pH was significant in all habitats and its influence on community composition remained significant even when habitat type was accounted for. Indeed, the important role of pH on bacterial communities has been observed on a global scale (Lauber et al. 2009), even if the direct mechanism by which it regulates microbial communities composition and functionality remains largely unknown (Malard and Pearce 2018). Lauber et al. (2009) suggested an indirect pH effect on the availability of different cations fundamental for life. A strong correlation between the relative abundance of genes encoding several metabolic and transport pathways and pH increase has been documented, suggesting a possible greater metabolic activity of bacterial cells in higher-nutrient and alkaline conditions (Bahram et al. 2018).

The C/N ratio, an indicator of substrate quality, was also an important parameter in determining changes in microbial communities (Table 1). N availability is an important determining factor for soil life forms (Chen et al. 2014; Leff et al. 2015). In general, oligotrophic species (k-strategists) dominate under N-limiting conditions, such as polar regions, while under abundant N concentrations copiotrophic species (r-strategists), able to utilize more labile C sources, prevail (Fontaine, Mariotti and Abbadie 2003; Chen et al. 2014). We found a higher richness and abundance of Actinobacteria and Armatimonadetes, and a lower richness and abundance of Planctomycetes in BSC plots than in BG and VV ones, with a higher C/N ratio (Figure S1, Supporting Information). Among these phyla, Armatimonadetes, although poorly studied, are generally considered oligotrophic (Lee, Dunfield and Stott 2014) and Actinobacteria, usually associated with plant roots, have been found in many desert soils (Anandan, Dharumadurai and Manogaran 2016), while Planctomycetes are usually more abundant in bulk soils than in the rhizosphere (Derakshani, Lukow and Liesack 2001). A higher relative abundance of Bacteroidetes and partially Proteobacteria, generally copiotrophic, were found in VV plots in respect to BG and BSC ones, with a more balanced C/N ratio (Figure S3, Supporting Information and Fig. 2).

Across the three habitats, the phylum Chloroflexi showed a clear preference for the bare ground habitat (Fig. 4b; Figures S1 and S3, Supporting Information). This phylum mostly includes oligotrophic organisms, apparently with greater stress-tolerance and/or lesser competitive capabilities than most other phyla (Costello and Schmidt 2006). In particular, the dominant Ktedonobacteria class showed a higher richness and abundance in the BG plots (Fig. 1 and Figure S2, Supporting Information) and was the class level grouping most influenced by edaphic parameters (Table S2, Supporting Information). This group is known to include organisms well adapted to extremely oligotrophic conditions and has been found dominant in volcanic soils of the Atacama Desert (Lynch et al. 2012) and in cinder deposits of the Kilauea volcano in Hawaii (King and King 2014) as well. Additionally, it has been recently shown that organisms belonging to this class are optimal CO and H₂ oxidizers and therefore considered pioneer organisms, allowing to fix atmospheric gasses in nutrient limiting environments (Islam et al. 2019).

Members of Verrucomicrobia, although still poorly studied, have been recorded in soils from many different biomes, even in Antarctica (Bergmann et al. 2011), and reported to be more abundant in plant rhizospheres compared to bulk soils in temperate environments (Jesus et al. 2010; Rocha, da van Elsas and van Overbeek 2010). Their close relationship with plants is also indicated by the fact that plant extracts must be added to culturing media in order to successfully isolate members of the Spartobacteria class (Sangwan et al. 2004). Although the richness of Verrucomicrobia did not differ significantly among the habitats studied, members of this phylum, including Spartobacteria, were significantly more abundant in BG plots (Figure S2, Supporting Information), possibly due to their greater tolerance to low nutrient conditions. In fact, this phylum includes many slow-growing organisms that are highly sensitive to changes in soil properties and, therefore, good indicators for changes in chemical factors linked to fertility (Navarrete et al. 2015).

Members of the phylum Cyanobacteria, despite the lack of a clear abundance pattern, had a higher richness in BG plots (Figure S1, Supporting Information). Filamentous cyanobacteria, in particular, have proved to be key organisms in early stages of soil development (Budel et al. 2016), due to their involvement in N fixation, moisture retention, soil surface stabilization and accumulation of organic matter in nutrient-limited environments, such as Antarctic ice-free regions (Cary et al. 2010). For this reason, we propose that BSC samples, characterized by abundant mosses and lichens, had a lower diversity of Cyanobacteria than BG plots because the latter reflects the earlier stages of microbial colonization.

CONCLUSIONS

This study provides a detailed picture of the landscape-level compositional dynamics of soil bacterial communities in Western Greenland, where information about soil microbiota is limited. It represents, to our knowledge, one of the first metabarcoding assessments of Arctic soil bacterial communities underlying different vegetation types. Bacterial richness did not correlate with increasing vegetation complexity, but there is evidence suggestive of a possible local loss of species connected to the expansion of shrubs at the expense of other soil habitats. Additionally, we found that community composition was strongly differentiated between the habitats and was strongly shaped by the vegetation composition in well vegetated plots.

The present report provides a status of the bacterial community composition, serving as a baseline for long-term monitoring, close to the vegetation transect ‘NERO line’ established to monitor future changes in the species composition of the plant communities.

ACKNOWLEDGEMENTS

Authors wish to acknowledge Hannele Savela (Thule Institute, University of Oulu, Finland), INTERACT TA Coordinator, for her suggestions and assistance in Arctic expedition, Claudia Pacelli (University of Tuscia, Italy) for her participation in the sampling activity, Christian Bay (Aarhus University, Denmark) for indications on sampling sites and Dr. Helga Szalontai, Mihály Jánószky (Eszterházy Károly University, Hungary) for the chemical analyses of the soil samples, and James T. Weedon for text language editing and his suggestions in the revision of the manuscript. József Geml acknowledges support from the MTA-EKE Lendület programme (no. 96049) of the Támogatott Kutatócsoportok Irodája. Fabiana Canini wishes to thank Naturalis Biodiversity Center for hosting her.

AUTHOR CONTRIBUTIONS

FC, LZ and JG planned and designed the experiment. CC collected the samples. FC performed DNA extraction. FC and JG performed data processing and analyses. FC, LZ and JG wrote the paper with inputs from SO and CC. FD actively participated in the revision of the manuscript.

FUNDING

This work was supported by INTERACT under the European Union H2020 Grant Agreement No.730938 (Project title: Effects of Climate change On Microbial Community of Soil in Greenland) and by the PRIN project 2015N8F555: Responses of alpine sensible ecosystems to climate changes (ReSaCC).

Conflicts of Interest

None declared.