Qiangyong Glacier is located on the north side of the Himalayas, with a length of 4.6 km, a maximum width of 2.8 km, and an area of 7.7 km² (Luo et al., 2003). It is a continental glacier; the snow line is at an altitude of 5,600 m, and the terminus is at an altitude of 5,000–5,100 m (Li et al., 2012). From 1976 to 2006, the glacier has retreated at an average rate of 4 m year^–1 (Yao et al., 2012). In July 2013, soil samples were collected from three sample sites (IS, SM, and BD) along two transects [line middle (LM) and line Top (LT)], on the east-facing side of the glacier movement residual moraine ridges (Figure 1). The three sample sites were as follows: the middle of the glacier termini and small Qiangyong lake upstream edge (IS), midway along the edge of the smaller lake (SM), and approximately two-thirds of the way along the edge of the larger lake (BD). At the time of sampling, sites IS_M(96a), SM_M(160a), and BD_M(259a) corresponded to a time of exposure after ice melting of 96, 160, and 259 years, respectively; sites IS_T(92a), SM_T(159a), and BD_T(258a) corresponded to a time of exposure after ice melting of 92, 159, and 258 years, respectively. The main plant located at IS and SM was Kobresia homilies, which then transitioned Potentilla fruticose at the BD sites. Bulk soil not in contact with the root system and located at 50–100 cm from each sampled abundant plant was collected from 3 to 5 cm after removing the top 1–2 cm of large sand grains. Soil samples near the plant root (less than 1 cm) were collected as rhizosphere soil from a similar depth as the bulk soils. Soil samples were stored in sterile sampling bags (Labplas, Canada) and transferred to the laboratory on ice. One part of each sample was stored at −80°C until DNA extraction was performed; the other part was used for physicochemical property analyses.

Soil factors were measured using the methods described in Guo et al. (2015) and will only be summarized here. Soil total C and total organic C (TOC) were measured in the solid state using a TOC analyzer (TOC-L, Shimadzu, Japan); soil total nitrogen (TN) was determined by elemental analyzer (vario MAX, Elementar, Germany). Soil nitrate (NO₃^–-N) and ammonium (NH₄⁺-N) were extracted with 1 M KCL and determined using an Automated Discrete Analyzer (AQ2, SEAL Analytical Inc., England). Soil water content (WC) was gravimetrically determined after drying at 105°C for 12 h.

Soil genomic DNA was extracted from 0.5 g of frozen soil using the FastDNA^® spin kit for soil (MP Biomedicals, Solon, OH, United States) following the manufacturer’s protocol. The quantity of the DNA was determined using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, United States). The primers 515F (5′-GTG YCA GCM GCC GCG GTA-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) were used for Illumina MiSeq sequencing at the Chengdu Institute of Biology, Chinese Academy of Sciences. PCR amplification and high-throughput sequencing method was according to Li et al. (2015).

The raw sequences were sorted based on the unique sample barcodes, and then quality control was performed for the sequences using QIIME pipeline (Caporaso et al., 2010). All sequence reads were trimmed and assigned to each sample based on their barcodes. The sequences with high quality (length > 150 bp, without ambiguous base “N,” and average base quality score > 30) were used for downstream analysis. Chimeric sequences were identified and removed using the Uchime algorithm (Edgar et al., 2011). After filtering and chimera removal, de novo operational taxonomic unit (OTU) picking was performed using uclust at 97% sequence identity, and subsequently, taxonomy was assigned to OTU based on the Greengenes database at a confidence level of 80% (version 13.8). All low-abundance Archaeal and Bacterial sequences that were only detected in one of the samples were culled. To avoid the influence of sequencing depth, rarefaction was performed with 9,000 sequences per sample for the diversity of microorganisms.

Community α-diversity was estimated using Shannon, Evenness, and Richness diversity indices. The β-diversity was calculated based on Bray–Curtis distances, and the differences in community structure between samples were visualized using principal coordinates analysis (PCoA) plot. Paired-samples t-tests were used for comparing data of bulk and rhizosphere soil with the SPSS 18.0 statistical software package. Adonis analysis was used to compare the difference between groups, and redundancy analysis (RDA) based on Bray–Curtis distances was used to determine the most significant environmental variables that might influence the prokaryotic community structure using the vegan and picante packages of the statistical platform R. Linear discriminant analysis (LDA) effect size (LEfSe) analyses between different sample groups utilized an online platform (Segata et al., 2011) and was applied to the OTU table to identify discriminant prokaryotic taxa.

The raw sequence reads generated in the present study are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the project ID PRJNA595717.

In total, the contents of TOC, TN, NH₄⁺-N, and NO₃^–-N in the rhizosphere were all significantly higher than in bulk soil (n = 18, p < 0.01). The TOC and TN contents showed an increasing trend along the glacier chronosequence, while NH₄⁺-N and NO₃^–-N contents showed a decreasing trend, at line middle and line top, respectively (Supplementary Figure 1).

In the phyla and genus level community composition, there were some clades detected in different abundance and succession trends between bulk and rhizosphere soil. In the phylum level, the dominant phyla (relative abundance > 5%) across all samples were Proteobacteria (21.95%), Actinobacteria (18.11), Acidobacteria (9.50%), Chloroflexi (7.9%), Crenarchaeota (6.81%), Bacteroidetes (6.27%), and Planctomycetes (5.69%), accounting for more than 75% of the total prokaryotic sequences. The relative abundance of the phyla Gemmatimonadetes, Verrucomicrobia, Euryarchaeota, and Firmicutes were > 1% (Figure 2A). The relative abundances of Proteobacteria and Bacteroidetes in rhizosphere soil was significantly higher than that in bulk soil; Crenarchaeota and Gemmatimonadetes were on the opposite (n = 36, p < 0.01). In bulk soil, the relative abundance of Acidobacteria showed an increasing trend along the chronosequence, while Actinobacteria displayed the opposite trend. A similar trend was also found in rhizosphere soil of the top line samples.

The dominant genera (average relative abundance > 0.5%) across all samples were Candidatus Nitrososphaera (5.7%), Arthrobacter (0.83%), Rubrobacter (0.78%), Modestobacter (0.55%), Flavobacterium (0.53%), Pseudomonas (0.53%), and Gemmata (0.52%). Seventeen other genera were found in the majority of samples with an average relative abundance of more than 0.18% (Figure 2B). Among them, the relative abundance of the predominant genus Candidatus Nitrososphaera in rhizosphere soil (3.69%) was significantly lower than that in bulk soil (7.89%) (n = 36, p < 0.001). The relative abundance of Modestobacter and Flavobacterium in rhizosphere soil was significantly higher than in bulk soil (n = 36, p < 0.001, p < 0.05). The relative abundances of Rubrobacter was on the opposite (n = 36, p < 0.05). Flavobacterium has an increased trend along the chronosequence. Modestobacter was decreased from site IS to SM and then increased from SM to BD in bulk soil, while in the rhizosphere they were on the opposite trend. Rubrobacter was relatively stable from site IS to SM, and then decreased from site SM to BD in bulk soil but increased in rhizosphere soil.

Prokaryotic community α-diversity indices showed mild succession in the glacier retreat chronosequence (Figure 3). In total, the averages of the three indices in the rhizosphere soil were slightly higher than those in bulk soils. In bulk soils, the indices were decreased from site IS to SM and then increased from SM to BD; the final value in BD was higher than that in IS in the line Top sample (Figures 3A,C,E). The opposite trend of succession was found in the rhizosphere soil in the line Top sample (Figures 3B,D,F). In line middle, the richness in bulk soil was relatively stable in bulk soil, but that has a drop in the rhizosphere soil, especially at site BD (Figures 3E,F).

Alpine plants affect the prokaryotic community β-diversity (Figure 4). Bulk and rhizosphere soils separated well from each other along the first axis (Figure 4A). This was in agreement with the Adonis analysis, which showed that the community composition of bulk and rhizosphere soils were significantly different (n = 18, p < 0.001). In bulk soil, each of the sampling sites (IS, SM, and BD) grouped separately (Figure 4B). In rhizosphere soil, the IS samples were close to SM along the positive direction of the second PCoA axis, but two subgroups were found in the BD samples, one of which (three samples all from the middle line) was much closer to the other two sampling sites (Figures 4C,E,F). According to Adonis analysis, the community compositions of IS, SM, and BD were significantly different in total, bulk, and rhizospheres soils, respectively (n = 12, 6, 6, p < 0.001). For each of the three sample locations (IS, SM, and BD), Adonis analysis showed that the bulk and rhizosphere soils within the same site were significantly different from each other (n = 6, p < 0.01) (Figures 4D–F). From IS to BD, the distance between the group centroids of the middle and top lines showed an increasing trend and no significant difference (n = 6, p = 0.576, 0.314, and 0.053, respectively) (Figures 4D–F). Together, these results suggested that the prokaryotic community composition was different at all three locations.

Across the glacier successional stages, different clades were detected in the bulk and rhizosphere soil fractions, which explained the statistically significant differences between their respective microbial communities (Figure 5). The numbers of discriminant clades were 24 and 11 from bulk and rhizosphere soils, respectively (Figure 5A). The number of discriminant clades increased in bulk soils from site IS to BD (6, 9, and 12), while the opposite trend was observed in rhizosphere soils (5, 4, and 3) (Figures 5E,F). In general, the phyla Proteobacteria, especially class Alphaproteobacteria, and Actinobacteria were significantly enriched in rhizosphere soil; meanwhile, Acidimicrobiia and Gemmatimonadetes were significantly enriched in bulk soil (Figure 5A). The number of discriminant clades in line top bulk and rhizosphere soils were both higher than in the respective line middle soils. In line middle, Alphaproteobacteria were significantly enriched in the rhizosphere soil, while Acidimicrobiia were enriched in bulk soil (Figure 5B). In line top, the phylum Proteobacteria, and specifically the class Alphaproteobacteria, was significantly enriched in rhizosphere soil. Meanwhile, Tenericute and Mollicutes were more enriched in bulk soil (Figure 5C). Among three sample sites in the total sample, site IS was enriched with Actinobacteria, SHA_109, and Thermoleophilia. Site SM was enriched with Betaproteobacteria and OP1.MSBL6. Site BD was enriched with Planctomycetia and Verrucomicrobia (Figure 5D).

RDA analysis was used to assess which factors affected the microbial community structure of the samples (Figure 6). In total, eight factors were detected: TOC, TN, NH₄⁺-N, NO₃^–-N, C/N ratio, WC, plants, and distance. For total soil samples, all factors except TOC had significant impact on the prokaryotic community structure (n = 36, p < 0.001) (Figure 6A). Among them, all the rhizosphere soils were significantly affected by plants. Distance affects both bulk and rhizosphere soil prokaryotic community structure; its effect on different samples sites was in the order of BD > SM > IS. The C/N affected the bulk soil microbial community structure more than that of the rhizosphere soil, while NH₄⁺-N, NO₃^–-N, and WC had the opposite trend (Figures 6A,D–F). In bulk soil, all factors except plants had significant effect on the prokaryotic community structure (n = 18, p < 0.001); the effect of C/N from different sample sites was in the order of SM > IS > BD (Figure 6B). In rhizosphere samples, NH₄⁺-N, NO₃^–-N, C/N ratio, WC, and distance were the main factors affecting the prokaryotic community structure (n = 18, p < 0.05); the effect of NH₄⁺-N, NO₃^–-N, and WC from different sample sites was in the similar order of IS > SM > BD (Figure 6C).

As expected, the special prokaryote community structures in the high-elevation Qiangyong Glacier foreland soil were revealed in this study. The relative abundance of seven phyla was above 5%; in a total of 11 phyla the relative abundance was above 1% (Figure 2A). The number of phyla was higher than that found in the frozen soil [five most abundant phyla (>2%)] from the glacier of the northwestern Himalayas of Jammu and Kashmir, India (Gupta et al., 2020); the two studies both were conducted in Himalayas. Furthermore, the three top abundant phyla detected in this research were also found in other glacier forelands of the Tibetan Plateau but with different order and abundance. For example, Proteobacteria (43%) and Actinobacteria (16%) were the predominant bacteria phyla in soils from Hailuogou Glacier retreat region. Similarly, they were also reported on the Muztag Ata Glacier chronosequence (Khan et al., 2020) and the frontier of Baishui Glacier No.1 (Sajjad et al., 2021). In this research, the Proteobacteria average relative abundance was 21.95% in total and comparable to the data from Gupta’s report (Gupta et al., 2020). The reason may be the effect of elevation, as in this research the glacier termini were at the elevation of 5,000–5,100 m; the height of the glacier is 4,700 m a.s.l. in Gupta’s report. The glaciers in other reports cited above were all below 4,400 m (Sun et al., 2016; Khan et al., 2020; Sajjad et al., 2021). These results suggested that as the elevation increased, the abundance of Proteobacteria decreased; at the same time, other phyla were enriched. As the microbial community reports on high-elevation glacier foreland were few, further research in this area is being looked forward to.

The relative abundance of Proteobacteria, Bacteroidetes, Modestobacter, and Flavobacterium in rhizosphere soil was significantly higher than that in bulk soil; some genera were in the opposite succession trend compared with bulk and rhizosphere soil (Figure 2). The LEfse detected differential clades infractions, which consistently explained the statistically significant difference between the bulk soil and rhizosphere soil prokaryotic communities (Figure 5). In general, the phyla Proteobacteria, especially class Alphaproteobacteria, and Actinobacteria were significantly enriched in rhizosphere soil; meanwhile, Acidimicrobiia, and Gemmatimonadetes were significantly enriched in bulk soil. At three sites, bulk and rhizosphere soil also had enriched clades. This result was different from reports found for low elevation, such as the finding that rhizobacterial communities were mainly composed of Acidobacteria and Proteobacteria, whereas bare soil was colonized by Acidobacteria and Clostridia in the foreland of Weisskugel Glacier (2,400 m a.s.l.) (Ciccazzo et al., 2014).

Meanwhile, in this research, the detected nutrient contents in rhizosphere soil were significantly higher than in bulk soil (Supplementary Figure 1). The α-diversity in rhizosphere soil was higher than in bulk soil (Figure 3); PCoA and Adonis analysis showed the bulk and rhizosphere soil microbial community was significantly different along the chronosequence (Figure 4). Based on the RDA analysis, plants also showed a significant effect on microbial community structures in Qiangyong Glacier foreland (Figure 6). These results suggest that the plants altered the soil C and N accumulation by shifting the soil prokaryotic community structure. This consisted of the most enriched phyla in rhizosphere soil (phyla Proteobacteria, especially class Alphaproteobacteria, and Actinobacteria) (Figure 5). Proteobacteria were detected as autotrophic microbial communities in Qiangyong Glacier-originated water (Kong et al., 2019). They were also supposed to have the ability for N fixation at a glacial foreland on Anvers Island (Strauss et al., 2012). Actinobacteria and Betaproteobacteria have the potential to photosynthesize; Betaproteobacteria and Alphaproteobacteria have the potential to perform nitrification and denitrification of microbes from supraglacial cryoconite of polar regions (Cameron et al., 2012). Furthermore, the dominance of Actinobacteria in Antarctic soils has been linked to specific trace elements (magnesium, calcium, and potassium) and salts in the associated glacier forelands (Bajerski and Wagner, 2013). Actinobacteria are adapted to oligotrophic environments where their hyphae allow them to restore nutrients and moisture through pores in the soil (Arocha-Garza et al., 2017; Zhang et al., 2019). The rhizosphere soil microbial development significantly affected soil organic C and total N accumulation in the Hailuogou Glacier forefield (Li et al., 2020). The results above suggest that the emergence of plants had a selective effect on microorganisms and promoted the accumulation of C and N in this high-elevation glacier retreat in rhizosphere soils.

The concentration of key nutrients related to soil fertility changed from site IS to BD along the Qiangyong chronosequence (Supplementary Figure 1). All detected environment parameters had effects on the high-elevation prokaryotic community composition succession (Figure 6). This result means that the distance, plant, and soil physicochemical properties together determine microbial community structures in the Qiangyong Glacier foreland. This result was partly similar to the reports of others. Microbial community structure is strongly conditioned by the successional stage, deglaciation time, water content, plants, spruce leachate, and the C and N contents of the foreland soil (Tscherko et al., 2004; Noll and Wellinger, 2008; Göransson et al., 2011; Knelman et al., 2012, 2018; Górniak et al., 2017; Kim et al., 2017; Bai et al., 2020). Most of the above reports did not include rhizosphere soil, and these studies were on elevations of < 1,000 m (Knelman et al., 2012, 2018; Kim et al., 2017), 1,780 m (Górniak et al., 2017), 1,920–2,054 m (Göransson et al., 2011), 1,950–2,050 m (Noll and Wellinger, 2008), 2,280–2,450 m (Tscherko et al., 2004), and 2,951 m at glacier termini (Bai et al., 2020). The effects of the high-elevation environmental parameters on the glacier foreland prokaryotic community were seldom found in field research. In this research, the effect of C/N on the bulk soil microbial community was greater than that in rhizosphere soil, while NH₄⁺-N, NO₃^–-N, and WC had the opposite trend. The above detected different physicochemical properties in each research, and most of them were in a relatively low elevation glacier retreat area. These results suggest that the distinct habitat was not ignorable. All the above results showed that plants, habitat, and glacier retreat chronosequence collectively control prokaryotic community composition and succession.