The physiochemical features of the glacier samples which are depth, temperature, pH, EC, Sulphate, nitrate, phosphorous, Potassium, Calcium–Magnesium, and pH ice and pH snow respectively, have been taken to define the physiochemical parameters, which were represented in the following Table 1.

From the Table 1, it was noted that temperature in winter was lowest in Biafo (− 20 °C), whereas it was maximum − 4 °C in the Kamri top. The pH of Burzil top was maximum (8.23), whereas the pH of Biafo was minimum (6.12). Electrical conductivity which was indicator of redox potential and metabolic activity too was highest in Biafo and Panmah. The nutrients concentration found highest in the Biafo where EC was also highest. Remaining sites showed less variability in the levels of nutrients present. Nutrient concentration as well as EC was considered as an indicator for prospect high diversity of bacteria in the sites. High temperature was considered to be indicator for hosting unusual bacterial species which were not representative organisms of glaciers.

The table also reveals seasonal variations in glacier physicochemical properties, with summer samples having shallower depths, lower temperatures, and slightly acidic pH values, while winter samples have higher electrical conductivity, sulfur content, nitrates, phosphorus, potassium, and calcium-magnesium content.

Raw data was subjected to quality assessment after trimming and truncation of low-quality reads. Trimming parameters were combined with standard filtering parameters, the most important being the enforcement of a maximum of 2 expected errors per reads (Fig. S1). After excluding non-quality reads all samples from G1 to G8 had more than 80,000 reads. Demultiplexed sequence was applied on the sequence forward reads as well as reverse reads and it observed that the minimum sequence count for forward read was 238nts, in forward reads its median was 238nt, and in reverse reads it was 250nts respectively. Per sample sequence counts of 10 samples showed that a higher number of sequence counts was observed in reverse reads sample (Fig. 2).

a shows the data obtained after demultiplexing sequences and parametric Box plot applied on the sequencing depth (parameter) on the samples, forward and reverse reads of the sequence obtained from the glacier samples, b The soil of the glaciers was tested and in forward reads the graph was generated using random sampling of 10,000 out of 695,238 sequences without replacement. The minimum sequence length identified during subsampling was 238 bases. Outliers are ignored for better results and symmetry. graph was generated using random sampling of 10,000 out of 695,238 sequences without replacement, c  OUT’s abundance was plotted in each sample, among all glacial samples G3 showed maximum OTU abundance and G7 showed minimum abundance.

The minimum sequence length identified during subsampling was 250 bases. The median was minimally influenced by outliers than the mean, and was generally the chosen central tendency indicator when the distribution was not symmetrical. Demultiplexing was the step involved in processing the information in order to know which sequences came from which samples after they had all are sequenced together. By applying the box plot on the obtained result on the samples taken from the 7 glaciers and one control, it was observed that the given samples had the regular linear increase in observed Phylogenetic diversity with increase in sequence depth, which means higher or larger the sequence depth more and greater the similarities (Fig. 2). Alpha diversity was analyzed to check the richness and evenness in different taxa in the seven glaciers. Operational taxonomic units were used to classify groups of closely related individuals against the sequence depth. Isotherm appeared when sequence depth was measured against each sample. Figure 2a is the isotherm which means there was a linear increase in observed OTUs with an increase in sequence depth, which means the higher or larger the sequence depth the greater the similarities. It is the regular linear increase in observed Phylogenetic diversity with the increase in sequence depth, which means the higher or larger the sequence depth, greater the similarities.

The alpha diversity of glacier samples labeled G1 to G8 was measured by Chao1, Shannon, Simpson, and Inverse Simpson metrics (Fig. 3). These metrics provided different perspectives on alpha diversity. Chao1 estimated total species richness, while the Shannon and Simpson indices provided information about the distribution of species and their relative abundances. The Inverse Simpson index also considered richness and evenness but placed more emphasis on rare species. By analyzing these metrics for these glacier samples (G1–G8), we got insights into the diversity, richness, evenness, and dominance patterns within each sample. This information gave us valuable for understanding for the ecological health and dynamics of glacier ecosystems.

Alpha diversity with reference to Chao1, Shannon, Simpson, and Inverse Simpson indexes. The Shannon index is in the range of 2.2 to 2.6. Generally, its value is between 1.5 and 3.5 in most ecological studies. The Shannon index rises as the community's richness as well as its evenness rises. This describes the index for a species' Beta diversity in a population since it is the combination of all the pairs collected from the various samples.

The y-axis values ranged from 140 to 200 represented the Chao1 values for each glacier sample. A higher Chao1 value indicated a higher estimated species richness in G3 sample. It indicates that G3 likely had a higher estimated species richness compared to other samples. Dots of G4, G5 and G6 were relatively close in height, it suggested that the estimated species richness was similar across these samples. On the other hand, there was a significant difference in dot heights for G1, G2, G7 and G8 indicating variations in species richness between these samples.

By observing the heights of the dots for Shannon values G5 suggested a more diverse community with a more even distribution of species with a Shannon value of 4.6. The dots of G3, G4, G6 and G7 were relatively close in height, which suggested that the diversity and evenness of species were similar across these samples. G1 was closer to G2 but there was a significant difference in values of G1 and G8, it indicated variations in species diversity and relative abundances between samples of those two glaciers. G8 and G3 showed a decreased value of the Shannon index which might be due to environmental stress.

The Simpson index measured the dominance or concentration of species in a community. It calculated the probability that two randomly selected individuals from the community belonged to the same species. A higher Simpson value indicated a lesser diverse community with a higher dominance of one or a few species. G5 and G7 had higher dots with a Simpson index greater than 0.92, which indicated that they had a higher dominance and concentration of species compared to G8 whose value was below 0.94 and G8 between 0.95 and 0.96. The Shannon index value of G1, G2, G4 and G6 were intermediate between 0.96 and 0.98 suggesting that the dominance and concentration of species were similar across those samples.

The Inverse Simpson index was the reciprocal of the Simpson index. It gave more weight to rare species, emphasizing their contribution to diversity. It quantified the effective number of equally abundant species. The y-axis values ranged from 20 to 80 represented the Inverse Simpson diversity values for each glacier sample. G5 and G7 had higher index values between 60 and 80 as compared to G3 and G8 whose values were below 20, which indicated that G5 and G7 had a higher richness and evenness of species, with a greater contribution from rare species. The dots of G1, G2, G4 and G6 were relatively close in height, which suggested that the richness and evenness of species, considering rare species, were similar across those samples. On the other hand, there was a significant difference in dot heights of G5 and G8, thus indicating variations in diversity, with more emphasis on rare species, between those samples. Overall G3 and G5 showed increased values of alpha diversity matrices thus showing an increase in diversity, richness, evenness, and dominance patterns. So overall highest abundance and diversity was seen in Siachin in summer and Kamri in winter.

Phyla Proteobacteria displayed higher OTUs accompanied by Actinobacteria, Firmicutes, and then Bacteriodetes in the heat map correlated with the relative percentage of each OTU in the samples. The Proteobacteria had the highest diversity in all the selected glaciers (Fig. 4a, b). Relative abundance in the percentage of total sequences at abundance of 0.01% presented distinct dominant phyla and specie in all samples. Distribution of OTUs were the taxa ’s relative abundance in the percentage of total sequences and the figure showed taxa with an abundance of 0.01 to 0.04%.

a, b  Heat maps showing the distribution of taxa at phylum, genus and species level. UC stands for uncultured OTUs.  NA stands for those OTUs which have not been classified yet. Values are the taxa's relative abundance in the percentage of total sequences and the figure shows taxa with an abundance of > 0.01 percent. All values are rounded to one digit; thus, the abundance in one sample of a taxon with a value of 0.0 lies between 0.00 and 0.04 percent. Different algorithms used gave slightly different results

Beta diversity analysis was used to measure and understand the variation in species composition between different samples within an ecosystem. Cumulative abundance of species along a gradient or ordination axis was estimated by CAP plots (Fig. 5a) which allowed us to visualize patterns of species distribution. In the above Fig. 5a CAP plot with and without aesthetics added was shown which referred to comparing two versions of a Cumulative Abundance Profile plot, one with basic representation and another with visual enhancements to better visualize and interpret patterns of species composition among different samples. The percentages on the x and y-axis indicated the proportion of explained variation for each axis.

a AP plot to visualize species distribution without aesthetics and with aesthetics added. b CA1 and CA2 with 44.9 coordinates at x and 36.0 at y-axis showed that control sample is lying differently as compared to all glacial sites. c Unifrac and unwunifrac analysis showed that with the distribution of samples as two subjects based on physichochemical parameters similarities and similarities in OTUs, both subjects again ended up in sharing same clusters

In Fig. 5b G1 was located farther to both axis. G2, G5, G1 and G7 were located in clusters indicated more similarity in species composition. Similarly, G4, G3 and G6 followed the same pattern. The greater distance of G1 from other samples showed its species dissimilarity from other glacier samples. In Fig. 5c samples were distributed in two subjects based on similarities in alpha diversity indices and shared OTUs. Both subjects were shared by all clusters which showed that variations among samples were not significant.

DCA was often applied as an initial step before choosing the appropriate ordination method. It was used to reveal the underlying structure in multivariate species abundance data (Fig. 6a). In the graph layout, the x-axis was labeled as DCA1 (51.6%) and the y-axis was labeled as DCA2 (23.4%). The graphs showed that dots were dispersed across the plot, which suggested a gradient in species composition and environmental factors. The G1 and G6 were present close along the y-axis while G7 and G2 were present close along the x-axis. It indicated similarity in their species composition and low scores on that axis. G5 was located farther away from the origin having high scores and showed dissimilarity from other samples composition. It showed that the species composition of kamri was less similar to other glacier samples. CCA plot was used to analyze and model the relationship between species abundance data and environmental variables (Fig. 6b). Different percentages on the x and y-axis indicated the proportion of variance for each axis. G1 being control clearly showed significant variation as compared to all other samples. Glacial samples clustered in two groups with G7, G6 and G2 being in one cluster and G3, G4, G5 and G8 in other cluster.

a DCA plot at DCA1 with 51.6% coordinates and DCA2 with 23.4% coordinates, b CCA plot with clustering of glacial samples into two distinct groupsc NMDS plot at NMDS1 and NMDS2 axis and Beta cluster hclust in the form of cluster dendrogram

Interpreting a NMDS plot involved understanding the arrangement of samples (dots) in the reduced-dimensional space represented by NMSD1 and NMSD2 axes (Fig. 6c). The positions of the dots provided insights into the similarity or dissimilarity of glacier species composition. In the above graph G1, G4 and G6 dots were closer to each other on the NMDS plot, which showed similarity in their species composition. G3 and G8 were also closer to each other. Similarly, G5, G2, G7 also followed the same pattern and depicted similar environmental conditions. Hierarchical cluster was used to identify patterns of community similarity and dissimilarity, which could have ecological implications and helped in the interpretation of species distribution across different environments or conditions. In the above Fig. 6c, the cluster dendrogram showed five clusters from top to bottom. There was a close similarity between G1 and G4. Glacier samples G3 and G8 had less distance and more similarity index. The same pattern was followed by G2 and G7. By moving from bottom to up more species added to form a new cluster i.e., G6, G1 and G4 showed similar species. Similarly, G2, G5 and G7 joined to form a cluster. At the top of the dendrogram all samples combined to form a single cluster. The reading scale showed that as more clusters joined, distance of vertical lines increased depicting less similarity index between species. For our understanding, we did cut the dendrogram at 0.10 reading to consider only those clusters that were more closely related.

PCoA on the basis of Bray–Curtis dissimilarity was performed by using the Bray–Curtis dissimilarity metric to quantify how different the species compositions were between pairs of samples. Then, we applied PCoA to these dissimilarity values to create a lower-dimensional representation of samples that preserved their pairwise dissimilarity relationships.

In Fig. 7a graph was plotted at axis1 with 56.3% coordinates and at axis 2 with 33.8% coordinates. The different axes of PCoA accounted for different variations. The result was in the form of two clusters based on phylogenetic distance metrics. PCoA focused on distance or dissimilarity matrices rather than raw data. The result was in the form of two clusters based upon phylogenetic distance metrics. G6, G1 and G4 were present in the form of clusters but at a distance from other glacier samples thus indicating greater pairwise distances between those samples. The G3, G7 and G8 were also clustered together indicating these groups of samples were more closely similar to each other than others. G2 and G5 were clustered together too. The clustering of dots on the graph also reflected species turnover which measured the change in species composition between different sites or conditions. Statistical tests, PERMANOVA helped to determine that observed clustering was statistically significant. RDA: CCA plot was a graphical representation to visualize the relationships between glacier samples and environmental variables. The above plots combined the strengths of both RDA and CCA to explore how environmental factors influenced the composition of species within samples. In the given figure factor datasets i.e., temperature, depth, Electrical conductivity, Nitrates, Phosphates, Sulphates, Ca–Mg, Potassium and pH were used to evaluate their effect on glacier samples (Fig. 7b). By looking at the graph it was observed that other than G3 all other sites OUT’s were directly impacted and influenced by the relevant physichochemical parameters suggesting strong redox potential and active metabolic activity at respective sites.

a PCoA plot on the basis of Bray–Curtis dissimilarity at axis 1 with 56.3 % coordinates and at axis 2 with 33.8% coordinates, b RDA:CCA plot at PC1 with 57.0% coordinates and at PC2 with 22.0% coordinates