Two mining complexes from Urvan and Progress were based on different soil types - Phaeozem and Umbric Gleyic, which determined variation in some soil chemical parameters, which we attribute to the region factor. Differences between disturbed and benchmark soils, which we classify as disturbance factor, are also reflected in some chemical parameters (Table 1). ANOVA with post-hoc Tukey HSD test showed that the region factor significantly affects soil temperature, pH, TC, potassium, ammonium and nitrates. The disturbance factor (between the benchmark and disturbed soil samples) affected OC and pH. Soil pH in the Urvan region was alkaline (7.4-8) and did not show a significant difference between quarry (UQ1 and UQ2) and benchmark (UB) sites (t = 1.09, p = 0.33) (Suppl. material 2, Table S2). Soil pH in the Progress region was slightly alkaline, ranging between 6.9 and 7.6, with quarry sites (PQ1-PQ3) being more alkaline (7.2-7.6) than benchmark PB soil (6.9-7.1) (t = 4.52, p < 0.01). OC quantities ranged from low to very low and had significant differences between quarry (0.2-0.8%) and benchmark (1.8-2.7%) sites for both regions (t = -14.60, p < 0.01). Some factors had significant correlation between each other (Table S3): phosphorus and ammonium (R² = 0.73, p < 0.05), ammonium and nitrates (R² = -0.53, p < 0.05), pH and nitrates (R² = -0.61, p < 0.05), phosphorus and nitrates (R² = -0.52, p < 0.05) and phosphorus and potassium (R² = 0.48, p < 0.05). Phosphorus content variation could not be attributed to the region or disturbance factor. Ammonium content was higher than nitrates in all sites, except PB (benchmark Agrisol) and PQ1 (2-year self-growing quarry with fresh dumps of soil and rock mixture), which is very close to PB. Both PB and PQ1 samples stood out from the rest as they had the smallest amount of phosphorus, ammonium and the maximum of nitrates. Interestingly, PB soil samples taken before and after the rain did not differ significantly in any parameters, except nitrates, which increased twofold after the rain. Samples from sites PQ2 and PQ3 from the same quarry bottom had the highest potassium content amongst all samples. To conclude, while pH, TC, ammonium, nitrates and potassium values demonstrated region specificity, OC values were associated with the disturbance factor, being higher in benchmark soils compared to primary soils of quarries in both regions.

A total of 84 libraries (four replicates for each of 21 soil cuts) of 16S rRNA gene were sequenced, resulting in 1768209 reads, which split into 10976 amplicon sequence variants (ASVs) or phylotypes. Minimum reads count per library was 7397, maximum 36229, mean - 21050. Reads not assigned on the phylum level (0.09% of the total count) were deleted from the dataset. From the remaining reads, 96.77% were attributed to class, 87.29% - to order, 70.03% - to family, 40.02% - to genus and 2.11% - to species. Microbiomes of different regions shared more than half (65.3%) of the total read count (Table 2). Microbiomes of benchmark and quarry sites in Progress had more common reads than those from Urvan (85.7% and 55.9%, respectively). The datasets generated and analysed for this study are available under accession numbers SAMN22858222-SAMN22858271 with links to BioProject accession number PRJNA777426 in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/).

Alpha diversity indexes, calculated for individual sites, had little variation, but some significant differences were observed. For both regions, we detected significant differences between benchmark and disturbed samples for the inverted Simpson index, which represents the probability that two randomly-selected sequences belong to different phylotypes (Fig. 2). In both cases, it was higher for quarry sites than for the benchmark. In Urvan, all indexes for separate samples from all sites differed significantly between each other with no relation to disturbance factor, with UQ2 site being the most diverse and UQ1 – the least (Suppl. material 1, Fig. S3). Apart from that, the dispersion of most alpha indexes of samples from Urvan was higher than in Progress. Indexes of alpha diversity allow us to estimate microbiome variation within samples and, in our case, we can assume that microbiome variation from Progress is more consistent across different sites and biological replicates, while from Urvan, it is more diversified.

a

b

Figure 2.

Alpha diversity in four indexes – Observed, Faith (PD), Shannon, Inverted Simpson for quarry/benchmark sites in different regions.

a: Progress region  
b: Urvan region  

Unlike alpha-diversity, beta-diversity revealed differences between samples in a more defined manner. PERMANOVA showed that microbiomes of quarry and benchmark samples differ significantly for both regions - with R² = 0.26 (p-value = 0.001) for Urvan and R² = 0.11 (p-value = 0.002) for Progress. Its values also show that disturbance is a greater factor of explained variability in microbiomes for the soils in Urvan than in Progress. Data visualised by NMDS matches with PERMANOVA, as we can see three distinct groups:

1. Urvan quarry samples UQ1 and UQ2,
2. Urvan benchmark UB and
3. All Progress samples PB, PQ1-PQ3 (Fig. 3).

Figure 3.  

Beta-diversity of soil sites. Biological and technical replicates of each site are surrounded by ellipses. Urvan: UQ1 and UQ2 – Quarry, UB – Benchmark. Progress: PQ1, PQ2 and PQ3 – Quarry, PB – Benchmark.

For the Urvan region, samples from both quarry pits group closer together, while benchmark samples are separated from them. On the other hand, for the Progress region, separation of microbiomes of benchmark and disturbed soils is less apparent.

Apart from differences between sites, we also investigated differences between biological replicates within one site. Analysis of multivariate homogeneity of group dispersions showed significantly higher distance to centroids values between replicates at Urvan sites (ANOVA p-value < 0.001) than for Progress sites (Suppl. material 1, Fig. S4). Thus, despite similar levels of differences in ecological microniches (variance in vegetation, climate and insolation), dispersion between biological replicates in Urvan was higher than in Progress.

The most abundant phyla across all samples were typical of soil microbiomes - Actinobacteriota, Acidobacteriota, Alpha- and Gamma- proteobacteria, Bacteroidota, Crenarchaeota, Firmicutes, Verrucomicrobiota, Planctomycetota, Chloroflexi, Myxococcota and Gemmatimonadota (Fig. 4). According to the heatmap, the relative abundance of phyla reflects differences between samples like beta-diversity: benchmark (UB) and quarry (UQ1, UQ2) samples from Urvan site demonstrate the difference in quantities of the phyla Bacteroidota, Crenarchaeota, Firmicutes, RCP2-54, Patescibacteria and Entotheonellaeota, while for Progress samples, there is no evident difference in major phyla composition between benchmark (PB) and quarry samples (PQ1-PQ3). As for location-specific phyla, Acidobacteriota, Planctomycetota and Chloroflexi show higher relative abundance in Urvan, while Verrucomicrobiota - in Progress.

Figure 4.  

Heat map of the phyla relative abundance across soil sites. Orange stands for the highest and blue - for the lowest values. Urvan: UQ1, UQ2 – Quarry, UB – Benchmark. Progress: PQ1-PQ3 – Quarry, PB – Benchmark.

On the family level, we still see that major groups are present in all samples, but their abundance differs between samples (Suppl. material 1, Fig. S5). Top taxa are Nitrososphaeraceae and Planococcaceae (higher values at benchmark sites), Chitinophagaceae (higher values at quarry sites), Pyrinomonadaceae, Sphingomonadaceae and Chthoniobacteraceae. Apart from these, Urvan sites have variation in the content of the following families between benchmark and quarry samples: Pseudonocardiaceae, Micromonosporaceae, Beijerinckiaceae, Bryobacteraceae and Comamonadaceae are more prevalent in the quarries UQ1 and UQ2, while Propionibacteriaceae - in the benchmark UB. For Progress, distribution of families does not seem to be linked to quarry/benchmark distinction, but rather to different quarry sites.

Statistically significant differences at the genus level between microbiomes from different sites were assessed by DESeq2 analysis, which estimates dependence between the mean read count (baseMean) and the variance of a phylotype between a set of samples. The outcome is expressed in a log2FoldChange value, which indicates how much the phylotype abundance has shifted between the compared samples. The higher the modulo value of log2FoldChange, the higher is the shift. Positive or negative values of log2FoldChange indicate the direction in which the shift occurs. We applied DESeq2 to detect shifts between regions (Fig. 5a), quarry and benchmark in Urvan (Fig. 5b), quarry and benchmark in Progress (Fig. 5c).

a

b

c

Figure 5.

Plots for DESeq analysis results. Dots represent phylotypes (ASVs), on the Y-axis is their baseMean and, on the X-axis, log2FoldChange value. The further the dot is from zero, the stronger the shift between compared groups, with negative values meaning more of the certain ASV in one group and positive - in the other.

a: Log2FoldChange values between sites from different regions, Progress on the left and Urvan on the right.  
b: Log2FoldChange values between quarry and benchmark sites in Urvan, quarry on the left and benchmark on the right.  
c: Log2FoldChange values between quarry and benchmark sites in Progress, quarry on the left and benchmark on the right.  

Only phylotypes with baseMean equal to 10 reads or more were left in the analysis, with log2FoldChange adjusted p-value < 0.05. With this cutoff region, comparisons revealed that there are 45 phylotypes, which are more abundant in Urvan and 164 - in Progress (Fig. 5a, Suppl. material 2, Table S4). The highest modulo values of log2FoldChange are detected for the phylotypes with baseMean < 100 reads, meaning that the highest differences are detected in minor phylotypes. On the contrary, phylotypes with baseMean values exceeding 100 reads have lower log2FoldChange values, meaning they are present in both regions, but most of them are prevalent in Progress. No apparent phylum tends to be more characteristic of any region. Both regions have different prevalent phylotypes from Acidobacteriota, Actinobacteriota, Bacteroidota, Crenarchaeota, Gemmatimonadota, Proteobacteria and others.

Comparisons between benchmark/quarry samples show different phylotype distributions in the two regions. In Urvan, there are 37 phylotypes which are more prevalent in the benchmark samples and 106 - in the quarry (Fig. 5b, Suppl. material 2, Table S5). The most abundant phylotypes in the Urvan benchmark site belong to Acidobacteriota (5), Actinobacteriota (6), Firmicutes (11) and Crenarchaeota (7). Phylotypes from quarry belonged to Acidobacteriota (24), Actinobacteriota (31), Bacteroidota (13), Proteobacteria (19) and Crenarchaeota (6). In Progress, four phylotypes were more abundant in benchmark soil microbiomes, all of which belong to Crenarchaeota; while 33 phylotypes were detected as more abundant in the quarry microbiomes, most of them belonging to Acidobacteriota (3), Actinobacteriota (10), Bacteroidota (8) and Proteobacteria (7) (Fig. 5c, Suppl. material 2, Table S6).

Several trends could be highlighted from the DESeq2 analysis data. Shifts of phylotypes quantities are the most contrasting between regions. Within regions, contrast between quarry and benchmark is more pronounced in Urvan than in Progress. Quarry microbiomes of both Urvan and Progress regions have a larger proportion of minor phylotypes in comparison to benchmark microbiomes. In all comparisons, phylotypes with higher baseMean values had lower modulo values of log2FoldChange, while phylotypes with low baseMeans have higher modulo values of log2FoldChange.

Canonical Correlation Analysis (CCA) was used to link beta-diversity of microbiomes from all sites with soil chemical properties. Biological replicates of microbiomes from different sites allowed us to build a significant model (ANOVA p-value = 0.002). For this analysis, temperature data were scaled to the deviation from the day’s mean to level the differences between sampling days. The most significant factors were pH (p-value = 0.003), OC (p-value = 0.009), ammonium (p-value = 0.006) and scaled temperature (p-value = 0.032) (Table S7). According to the VIF test, these factors were not multicollinear, meaning their values cannot be linearly predicted from each other and their effect on sample microbiomes in this model was independent (Suppl. material 2, Table S7). Data were visualised on the CCA plot, where arrows point the direction of the factors' influence and dots represent microbiomes or individual phylotypes (Fig. 6). The further the dot is relative to the direction of the arrow, the more the factor explains the variation. If the dot is in the direction opposite to the arrow, then the value of the factor correlates negatively with the composition of the microbiome or the presence of the phylotype. In Urvan, microbiomes of quarry sites UQ1 and UQ2 are mostly affected by pH and ammonium, while those of benchmark UB site are mostly affected by OC. All microbiomes from the Progress site are influenced by nitrates and TC quantities. Progress samples show less dependence on chemical factors than Urvan because they are grouped closer to each other and are closer to the 0.0 point. On the other hand, Urvan samples are quite dispersed far from 0.0. In accordance with the beta-diversity, microbiomes from benchmark and quarry samples from Urvan (UQ1 and UQ2) exhibit more pronounced differences between each other, which can be linked to the influence of chemical parameters, than samples from Progress (PQ1-PQ3), which group closer to each other (Fig. 6a).

a

b

Figure 6.

Plots for CCA analysis, the further the dot is from the arrow, the more it is influenced by the factor. Dot - a microbiome of a single soil cut or a single phylotype. TC – total carbon, OC – organic carbon, temp – temperature scaled to day’s mean.

a: Relationship of the microbiome composition of individual sites to soil chemical parameters.  
b: Relationship of the top 100 abundant phylotypes (expressed in ASVs) to soil chemical parameters.  

The CCA plot on Fig. 6b shows how singular phylotypes in this dataset are affected by chemical factors. The majority of the top 100 abundant phylotypes are dispersed in the direction of OC and nitrates. These include phylotypes from different phyla, the most reactive ones belong to Bacillus, Nitrososphaeraceae, Microlunatus, Acidibacter, Xiphinematobacter, Nitrospira and Gaiellales. Some of these phylotypes match those which are statistically more prevalent in benchmark sites (Suppl. material 2, Table S5 and S6). In the opposite direction of OC, there are phylotypes, associated with the lack of OC. These are RB41, Ramlibacter, Flavisolibacter, Asanoa, Puia, Niastella and Briobacter. Some of these phylotypes are also detected as characteristic of quarry microbiomes by Log2FoldChange analysis. Other factors influence microbiome composition, which could be linked to some other differences between samples, not linked to disturbance or region factors. Such factors include ammonium and nitrates: they influence the microbiome composition in opposite directions, so that phylotypes reacting to the presence on ammonium (Pseudonocardia, Solirubrobacter, Acidibacter) and nitrates (Udaeobacter, Nitrososphaeraceae, Gaiellales) could be distinguished.

We also used CCA to show the connection between microscale (chemical parameters) and macroscale (region, disturbance) factors and their influence on microbiome composition. If we put on the CCA plot only phylotypes, significantly changing between regions (detected by DESeq, Suppl. material 2, Table S4), we see that phylotypes split into two groups: the ones from Urvan are on the side of ammonium and pH influence, while on the opposite side, there are phylotypes from the Progress region, which are more influenced by TC and nitrates. Both regions have characteristic phylotypes from Acidobacteriota, Actinobacteriota, Proteobacteria, Crenarchaeota and other phyla.

Fig. S6b (Suppl. material 1) demonstrates the CCA plot with phylotypes that were revealed as statistically different between quarry and benchmark samples for each region by DESeq analysis. Here, phylotypes are located on opposite sides of OC, which is consistent with the fact of correlation between OC quantities and soil disturbance. These phylotypes are mostly from Acidobacteriota, Proteobacteria and Actinobacteriota. Groups corresponding to the presence of OC (Crenarchaeota, Firmicutes) are divided into two: one between OC and ammonium, corresponding to phylotypes from the Urvan benchmark and the other between OC and nitrates, corresponding to the Progress benchmark. Notable, phylotypes with the same taxonomy (Nitrososphaeraceae, Vicinamibacteraceae, Gemmatimonadaceae) can be seen in the groups correlating with both the presence and the absence of OC.

Analysis of 16S rRNA gene libraries is often based on the negative binomial distribution, which has some limitations. These methods make a Type I error in assessing changes in the microbial community at high taxonomic levels (Lin and Peddada 2020, Nearing et al. 2022). The PhILR transformation offers an approach to overcome statistical artifacts of relative abundance of microbiota and analyse compositional data. It reveals “balances” which allows further investigation into the relating of phylogenetically close microorganisms to different factors. We applied this approach to detect changes in microbiota associated with the disturbance factor from each region. Nine significant balances for Urvan and seven significant balances for Progress were identified. It was shown that, for Progress, significant differences between benchmark and disturbed soils appear at low taxonomic levels (differences in individual phylotypes within genus or family) (Suppl. material 1, Fig. S7A). In contrast, for Urvan, several balances (n210, n733, n788) show differences at the phylum-class levels (Suppl. material 2, Fig. S7B). The response is especially diverse at different taxonomic levels within the Acidobacteriota phylum. Differences in the phylotype content within the phyla Firmicutes, Chloroflexi, and Verrucomicrobiota is also shown in the quarry sites. These results support previous observations: microbiomes of disturbed, but reclaimed quarry soils in Progress are much more similar to benchmark soil than microbiomes of unreclaimed quarry soils in Urvan.

This work is a continuation of the research on the microbiomes of soils from different climate zones, recovering from anthropological damage, primarily mining of parent rock material (Gladkov et al. 2019, Ivanova et al. 2020, Pershina et al. 2020, Zverev et al. 2020, Abakumov et al. 2021, Kimeklis et al. 2021). Here, we explore the area with a moderate continental climate, located in the northern foothills of the Caucasus mountains. The specificity of this study is that we collected samples from the same type of sandy-gravel quarries in close regions (approximately 50 km) with different benchmark soil types. The key difference between sampling sites was that quarry pits at one region (Progress) were reclaimed by backfilling using soil heaps from the overburden Agrisol, while the others (Urvan) were left to passively recover with spontaneous vegetation overgrowth on parent rock material. Thus, we were able to analyse different patterns of microbiome restoration on similar parent rock material in one climatic zone, but with different benchmark soil types and applied reclamation practices.

In our previous studies from the quarries of northern regions, we observed that primary soils of quarries are colonised by photosynthetic bacteria – Cyanobacteria and Chloroflexi (Gladkov et al. 2019, Kimeklis et al. 2021). These microorganisms form biofilms or soil crusts and can successfully colonise substrates deprived of organic carbon resources (Malard and Pearce 2018), but in this work, presence of these groups of microorganisms was minuscule. Perhaps, it can be explained by the overall lower humidity of the southern region and warm arid conditions during the period of sample collecting.

Factors that influence microbiome composition can be put into hierarchical categories, based on their complexity and scale of effect (Deakin et al. 2018). Large-scale differences, like regions, distance or type of agricultural practice, usually are considered the main factors of microbiome variation (O'Brien et al. 2016, Deakin et al. 2018, Shi et al. 2018). On the other hand, local overgrown vegetation, which is one of the bases of humic layer accumulation (Abakumov et al. 2020), usually creates spatial variations, which translate into microscale differences in microbiomes (Schreiter et al. 2014, Mitter et al. 2017). The same goes for comparisons between seasons and distance - while seasons create variation due to fast-changing factors, geographical differences explain more of microbial variation (Zhang et al. 2020, Wang et al. 2021). In our study, we detected differences between biological replications caused by microscale factors, like vegetation, insolation or water regime, but they did not overcome diversity created by large-scale factors - region and disturbance. There is evidence that the application of soil reclamation techniques shortens the recovery period and stabilises the microbial community (Hou et al. 2018). Here, we detected the same effect: application of backfilling in the quarries of the Progress region led to a significant reduction in microbiome dispersion and the difference between disturbed and benchmark sites. On the other hand, quarry microbiomes of the Urvan region demonstrated higher distinction from benchmark samples and higher dispersion of biological replicates. This effect can be explained by the fact that, in poor unreclaimed gravel heaps, microbiota has higher sensitivity to microscale spatial variation of nutrients introduced by plants than in primary soils mixed with Agrisol (Naylor et al. 2020, Ayangbenro and Babalola 2021).

Another effect that happens with soil disturbance is the adaptation of microorganisms to the new conditions. It was shown that, in treated soils, relevance of abundant microorganisms (bacteria and fungi) is reduced and the relevance of low-abundance microorganisms is increased (Bossolani et al. 2021), which can happen since conditions have become less favourable for major microbiota and more favourable for the growth of minor microbiota. We observed this effect in both regions: in comparison to benchmark soils, in disturbed soils, major phylotypes decrease their abundance, while higher quantities of minor phylotypes emerge. In Progress, microbiomes of disturbed soils still carried the same major phylotypes from Firmicutes and Crenarchaeota as in benchmark Agrisol, but their quantities were relatively lower. At the same time, disturbed soils carried many minor phylotypes, some of which from cellulose decomposing Chitinophagaceae and Cellulomonas. The presence of these taxa can be linked to the increased content of plant residues, which accumulate in quarry soils due to lack of crop harvesting (Kolton et al. 2013).

The content of the most measured chemical parameters, including OC, phosphorus and nitrates was low across all sampling sites, which is typical for the local soils (Gorobtsova et al. 2017, Gorobtsova et al. 2021). The acidity of the benchmark and quarry soils corresponds to their genetic features, which is explained by the chemical composition of the mineral waste (Gorobtsova et al. 2016, Gorobtsova et al. 2017, Gorobtsova et al. 2021). In Urvan region, the main difference between benchmark and quarry soils was reflected in carbon content: benchmark soil retained higher percentages of organic carbon, while primary soils of quarries showed high quantities of total carbon, enhanced by parent rock material. Soil cover in Progress reacts differently to the introduction of parent materials: initial Agrisol already has high quantities of carbon, it rises in the freshly reclaimed quarry bottom, but is reduced in older quarry bottoms. Organic matter content reduction in disturbed lands with reclamation was reported earlier (Liu et al. 2017).

Benchmark soil in Progress - Agrisol - was the only soil showing the prevalence of nitrates over ammonium, which is typical for agricultural soils (Cui and Song 2007). This feature is preserved in the mining pit neighbouring the field, which showed signs of recent reclamation with the mixture of rock dumps and soil heaps. In the soil of the older quarry, the balance of nitrates and ammonium shifts back to ammonium prevalence. It could be linked to soil acidic status, nitrate leaching or an introduction of plant residues since crops were no longer being harvested (Dejoux et al. 2000, Miller and Cramer 2005).

Using several biological replicates with varying chemical parameters allowed us to create a reliable model of factors influencing the microbial community. The key factor defining soil disturbance in both regions was OC, which revealed the same phylotypes from the two regions reacting to its content - Nitrososphaeraceae in benchmark soil and Nitrososphaeraceae, Azospirillaceae, Cellulomonadaceae, Vicinamibacteraceae, Nocardioidaceae and Chitinophagaceae in quarries. Members of Azospirillaceae were reported to be associated with plants and to be involved in carbon and nitrogen cycles (Sun et al. 2020). Vicinamibacteraceae from Acidobacteriota were reported to be flexible in their preferred carbon source (Navarrete et al. 2015). Cellulomonadaceae can degrade not only plant residues, but other carbon sources, like DNA and chitin (Stackebrandt and Schumann 2014). Nocardioidaceae are considered mostly chemoorganotrophs (Tóth and Borsodi 2014). Thus, with the lack of easy organic carbon in disturbed soil, microbiomes are becoming enriched by microbiota, flexible to available energy resources.

Traditionally, microorganisms are divided by their life-history strategy into fast-growing copiotrophs (or r-strategs) and slow-growing oligotrophs (K-strategs) (de Vries and Shade 2013). Based on their growth rates, usually gram-minus soil bacteria, like Proteobacteria, Bacteroidota and Gemmatimonadota are treated as copiotrophs, while gram-plus bacteria (Firmicutes, Actinobacteriota) as oligotrophs (Fierer et al. 2007, Zhang et al. 2020). However, this division is quite arbitrary and does not always follow taxonomic division (Ernebjerg and Kishony 2012, Ho et al. 2017, Song et al. 2017). For example, archaea from Nitrososphaera are reported to correlate with nitrate composition in soil (Zhalnina et al. 2014), but in our dataset, we found different Nitrososphaera phylotypes correlating with contents of nitrates, ammonia and organic carbon. Moreover, Ramlibacter representatives were described from poor nutrient desert environments (Heulin et al. 2003) and alongside this fact, we found phylotypes from Gammaproteobacteria - Ramlibacter and Ellin6067 - associated with the lack of OC. However, there were also other phylotypes attributed to Ramlibacter and Acidibacter, associated with the presence of OC. The same trend was detected in Acidobacteriota - phylotypes from Vicinamibacteraceae and Blastocatellaceae were detected in OC-rich soil samples and Vicinamibacteraceae, Bryobacter and RB41 in OC-deprived soils. Thus, while the prevalence of a certain phylum in the dataset can be linked to some microbiome-forming factors, our analysis showed once again that phyla are formed by phenotypically paradox lower taxons.