Soil microbiome of different-aged stages of self- restoration of ecosystems on the mining heaps of limestone

Soil microbiome plays an important role in soil forming process as well as soil functioning. This is especially relevant for initial stages of soil regeneration after strong anthropogenic impact (i.e., in quarrying complexes). The study of soil microbiome in such areas is crucial for both understanding the forces driving soil formation and optimization of the reclamation techniques. This study is aimed to investigate the soil microbiome of different ages of soil formation on the heaps of limestone mines in conditions of boreal climate of North-West of Russia. Soil microbiome analysis was performed using high-throughput sequencing. Soils of the investigated sites were predominantly Rhendzic Leptosols. The initial set (2 years) ofOTUs in the entire period of soil formation remains unchanged, accumulating more than 98% of the total microbiome in each time point. Analysisofbeta-diversity (weightedandunweighted)demonstrated a clear differentiation of microbiomes of all time points, which suggests that underlying taxonomic structure variations can be attributed to variation of the minor taxa set composition. The most significant differences occur in Proteobacteria, Chloroflexi, Acidobacteria, and Actinobacteria, butmajorityof thesephylahavedifferentsetsofOTUsdemonstrating opposite trends. Generally time positively correlate withmicrobiome of carbonate soils in all cases of benchmark and newly formed soils is familiar because of the uniform zonal environmental conditions in geochemically specific landscape of Izhora upland.


Introduction
Soil microbiomes play an important role in soil functioning and implementation of key ecosystems services (Sprocati et al. 2014; Abiraami et al. 2019; Gorzelak et al. 2020). Nowadays, the possibilities of new generation sequencing allow to receive huge amounts of data on taxonomical and functional diversity of soil microbiome . These data can be used for management of microbiological drivers of soil formation. Key soil forming processes are connected with activity of soil microbiome directly or indirectly. This is especially important for the very initial stages of soil regeneration after strong anthropogenic impact (Urbanová et al. 2011;Sprocati et al. 2014). That is why an accumulation of the data on soil microbiome in various postanthropogenic and technogenic environments is quite useful for creation of fact-based databases for environmental engineers and decision makers in the field of soil and vegetation cover restoration (Dmitrakova et al. 2018a(Dmitrakova et al. , 2018b. At the same time, we have now more raw and processed data on microbiome of zonal mature soils (Taş et al. 2018) than on microbial communities of abandoned or reclaimed lands, recently involved to open-cut mining (Urbanová et al. 2011). It is quite important to understand the structure of microbial communities on post-mine spoil heaps, both for understanding the forces driving soil formation and applications in the optimization of the reclamation techniques. Soil restoration in an initial stage is driven by a few groups of microorganisms, but not all of them are presented by cultivated forms, especially on very initial stages of pedogenesis. This is especially important for initial stages of abiogenicbiogenic interactions in soil and soil-like bodies, formed on the surfaces of recently active mines. New generation sequencing allows investigating completely new levels of biodiversity of anthropogenic environments than classical methods of laboratory cultivation of microorganisms ). Previously, this method was recognized as an informative and productive tool for investigation of soil microbiome on the mines of various bioclimatic zones. Thus, huge data were obtained and published for mining heaps of polar environments (Gladkov et Mokma et al. 2004). Essential parts of mining-affected areas in Russian North-West and in adjacent Estonia are presented by open mines of limestones, phosphorites, and shales. The technologies of soil reclamation are considered as outdated practices. They should be amended and updated by methods, based on increasing or management of microbiological activity of developing soils on various stages of the restoration process. Thus, the aim of this study was to conduct for the first time an investigation of soil microbiome of different ages of soil formation on the heaps of limestone mines in conditions of boreal climate of Russian-North West.

Study area
The study sites are located in the central part of Leningrad region (North West of Russia, Figure 1), on Izhorskaya upland (another name -Odrovician plateau). This is an upland which has been used for intensive agricultural practices in time of Sweden's presence in this area due to high fertility of soils with high percentages of carbonates -Rhendzic Leptosols (FAO 2014). Thus, this region was intensively anthropogenically and industrially affected. In the beginning of the nineteenth century, intensive mining of limestones was started with the aim of supplying Saint-Petersburg city with limestones for buildings and facing stone. One of the most ancient limestone quarries in Leningrad region is located in Pudost settlement and it was designated for limestone mining for Kazanskiy Cathedral construction and facing. Elizavetino quarry is the second largest and oldest mine in Izhorskaya upland. The age of this quarry is about 200-220 years, but soils of the heaps are younger, because of many replacement materials of heaps accumulated inside the quarry. There are numerous similar quarries, located in the west direction, till the border with Estonia  following herbal species: Fragaria vesca, Calamagrostis officinales, Vicia cracca, Tussilago farfara, Melilotus albus, and Gallium album. The tree species are presented mainly by alder (Alnus incana) and willow (Salix caprea) with some admixture of apple tree (Malus ideaus). At the stage of 70 years, the developed rare-stand forest is composed of willow, alder, poplar, rowan (Sorbus aucuparia), aspen (Populus tremula) with mixed and rich herb stand, and essential part of bryophyta. The benchmark soil is located under the zonal developed mature deciduous-coniferous (Picea abies, Tilia cordata, Acer platanoides) forest with prevalence of nemoral species (e.g., Asarum europaeum and numerous ephemeroids). This soil can be classified as Rhendzic Leptosol, formed on derivatives of limestones.

Laboratory analyses
Soil samples for top horizons were sampled into the plastic bags and frozen in the field; other sample types were taken from every soil horizon with the aim to characterize better soil profile development rate and vertical soil profile structure. These unfrozen samples were transported to the laboratory, grounded, and passed through 2 mm sieve.
The particle size distribution was determined using the Kaczynski sedimentation method with pyrophosphate peptization of microaggregates (Rastvorova 1983), the organic carbon content was determined by C-N-H analyzer (EA3028-HT, Research Park of Saint-Petersburg State University), and the pH of the aqueous extract was determined potentiometrically in a soil:water ratio of 1:2.5. The concentration of mobile phosphorus and potassium compounds was measured as per Kirsanov in TsINAO modification, exchangeable ammonium by TsINAO method (GOST 2011), and nitrates ionometrically (GOST 1985(GOST , 1986. DNA was extracted in 6 replicates for each site from 0.5 g of soil using the NucleoSpin ® Soil Kit (Macherey-Nagel GmbH & Co. KG, Germany) using a combination of SL1 + SX buffers. Homogenization of soil samples was performed on a Precellys 24 device (Bertin Technologies, France). The purity and amount of DNA in the preparation was assayed electrophoretically in 1% agarose with 0.5 × TAE buffer (DNA concentration in a sample averaged 50 ng/ml). Purified DNA was used for quantitative PCR (qPCR) and preparation of amplicon libraries as per the instructions of the sequencing protocol (Illumina, Inc., USA).
To estimate the quantity of bacteria in the soil samples, quantitative qPCR was performed with the following primer pairs for the 16S rDNA fragment -EUB338 (5′-ACT CCTACGGGAGGCAGCAG-3′) (Lane 1991) and EUB518 (5′-ATTACCGCGGCTGCTGG-3′) (Muyzer et al. 1993). qPCRmix-HS SYBR kit (Eurogen, Russia) was used to prepare the reaction mixture according to the manufacturer's instructions. The standards were the series of 10-fold dilutions of fragments of the 16S rRNA gene of E. coli. Each sample of the PCR mixture, including standards, was analyzed in triplicate. The measurements were carried out on a CFX96 thermocycler (Bio-Rad, Germany) according to the following protocol: 3 min at 95°C, 20 s at 95°C, 20 s at 50°C, and 20 s at 72°C (40 cycles). For replicates of both PCR and different DNA samples in a site, the mean values (M) and standard errors of the mean (±SEM) were calculated. After processing, the results were expressed as the number of ribosomal operons per 1 g soil. Significance of mean differences was assessed by Student test with p < 0.05.
In preparing libraries of 16S rRNA gene fragments for each soil DNA sample, PCR was performed with universal primers for the variable region V4, the F515 5′-GTGCCAG CMGCCGCGGTAA-3′ and R806 5′-GGACTACVSGG-GTATC TAAT-3′ and QIIME (Caporaso et al. 2010) software according to the following process. All sequences were trimmed, paired-end reads were assembled, and the quality of nucleotide sequences was checked. All nonbacterial and chimeric sequences were removed, and data were normalized. Sequences with more than 97% similarity were combined into operational taxonomic units (OTUs) using the de novo UCLUST-based algorithm. One sequence was selected from each OTU to compile a set of representative sequences. The representative sequences were classified using RDP naïve Bayesian rRNA Classifier program, and then aligned according to PyNast algorithm, using the SILVA database as a matrix.
For a comparative analysis of communities, the indices of alpha-and beta-diversity were calculated. Alpha-diversity was evaluated using species richness indices (the number of OTUs in the sample (Observed)), Faith's phylogenetic diversity index (PD), and the Shannon and Simpson indexes of diversity (Faith, 1992).
The significance of differences between microbiomes by alpha-diversity indices was evaluated using Mann-Whitney test. To determine beta-diversity, weighted and unweighted unifrac methods were used. The results were analyzed using multivariate statistics (analysis of the principal coordinates with R analysis). A Mantel test (Pearson correlation, 100 permutations) for combined replications was performed with QIIME to assess the correlation of bacterial community composition with concentration of basic soil macroelements (Mantel 1967). To assess the significance of differences between abundance of individual taxa in different soils, DeSEQ2 library (Love et al. 2014) was used.
The paper presents the mean values of indicators (M) and their standard errors (±SEM). The revealed differences were considered statistically significant at p < 0.05.
Ethical approval: The conducted research is not related to either human or animal use.

Soil development and chemical features
Soils, formed on the bottom of quarry, represent not a fast mode of initial soil development due to the fact that they are formed on the carbonate containing materials. In such cases, the speed of soil evolution is lesser than that on acid and weathered parent materials Data on soil chemical composition are provided in Table 1. These data illustrate that soils investigated are typical representatives of Calcaric soils in the initial or primitive branch of national soil classification (2004). These soils are alkaline, and the degree of alkalinity decreases with the age as well as carbonate content. This is quite typical for different-aged soils of differentaged ecosystems of Ordovician plateau (Blagovidov 1946;Gagarina et al. 1981;Gagarina 1996). Decreasing of carbonate content during leaching and weathering results in decreasing of skeletal fraction content and increasing of the fine earth percentages which well corresponds with data of Abakumov et al. (2005).

Soil microbiome
Data on soil microbiomes are presented on Figures 3-6. The distribution of copy number of 16S rRNA operons per 1 g of soils on the different age stages of soil formation ( Table 2) shows the increasing trend within the age row. This parameter indicated the bulk level of biologization of substrate from the initial stage of soil formation till the climax stage, e.g., benchmark ecosystems. These data are in good correspondence with metagenomics data, obtained previously (Urbanová et al. 2011).
The analysis of the microbiome dynamics in the course of soil formation has shown several important features. As expected, the total sum of bacteria is constantly increasing through 4 soil ages starting with very low 10 7 and ending with 10 9 at the "mature" soil benchmark site ( Table 2). This is in complete accordance with the soil carbon/nitrogen accumulation which is clear mark of the soil formation where bacteria are one of the most important players. But something is surprising: as it is evident from the Figure 3, the taxonomic composition (on the phylum level) remains almost the same through all 200-years process of the soil formation. Even the initial OTU set of the first age (2 years) remains unchanged all the history of the soil formation: the fraction of the total microbiome contained in the initial set of OTUs in all soil samples corresponding to 26, 70, and benchmark soil is about 98% in all cases. The alphadiversity indices describing richness, evenness, and phylogenetic diversity also vary in a very narrow range (Figure 4), and according to the statistical tests, have no significant difference between soil ages. But, behind a stable picture of the alpha-diversity, a real dynamic process can be detected based on the weighted and unweighted betadiversity analyses ( Figure 5). From the data presented, it is clear that mechanisms of the soil microbiome evolution in course of the soil formation include both changing the taxa set (unweighted unifrac) and their abundances (weighted unifrac). Interestingly, both ordinations demonstrate the close similarity between soil communities of 26 and 70 years samples and separate location of the initial and final points (2 years and benchmark). These data taken together demonstrate that the evolutionary rate of the soil microbial community is relatively high in the first years of the soil formation and has trends to a slowdown in the course of further evolution.
There are some important details about weighted and unweighted beta-diversity metrics. Based on the abovementioned finding that about 98% of the total communities of all soil samples of all soil ages belong to the initial set of OTUs detected in the 2-years old samples, one can propose that the dynamics of the taxonomic structure detected in the unweighted Unifrac could be attributed only to the minor prokaryotic taxa. Contrary to the unweighted metric, the difference between the soil microbial communities evident from the weighted unifrac ordination can be attributed to the taxa representing an essential part of the microbial community. This is evidenced by both clear separation of the microbial communities of different ages and the high level of the explained variabilities (axis 1 and axis 2) which are twice the unweighted values.
For revealing active taxa (which have positive or negative abundance changing), we use DeSEQ2 algorithm. By beta-diversity, we can see nonlinear changes in communities: two clusters (22 and 70 years soils) are, in fact, really close. So we perform DeSEQ2 analysis only for two pairs: 2-26 years samples and 26-200, respectively. According to DeSEQ2 comparison, most significant differences occur in Proteobacteria, Chloroflexi, Acidobacteria, and Actinobacteria. Those phyla contain significant amounts from all active communities. Meanwhile, many OTUs from      . Surprisingly, the initial set (2 years) of OTUs in the entire period of soil formation remains unchanged, accumulating more than 98% of the total microbiome in each time point. This can be explained only by assuming that at the initial stage of the soil formation the initial mineral material is inoculated with the whole microbiome set transferred from the adjacent benchmark soil sites probably with aerosols. Due to that, observed soil microbiome dynamics cannot be used as a model for the primary soil formation in the past when the "benchmark inoculum" was unavailable. On the other hand, that suggests that mineral material inoculation with benchmark soil microbiome can direct the soil formation in the required direction and could be used in the real practice of the soil restoration. Moreover, it is necessary to take into account that the landscape of Izhora upland is intrazonal -Rhendzic Leptosol, so the microbiome of benchmark soil could be familiar with very initial soils due to pedoenvironmental specificity. 4. Despite the initial OTUs set conservation during the whole chronoseries, analysis of beta-diversity (weighted and unweighted) demonstrated a clear differentiation of microbiomes of all time points. This suggests that underlying taxonomic structure variations can be attributed to variation of the minor taxa set composition (explains unweighted beta-diversity) and variation in abundances of major taxa (explains weighted betadiversity). 5. Analysis of the statistically significant changes in taxa abundances shows that the most significant differences occur in Proteobacteria, Chloroflexi, Acidobacteria, and Actinobacteria, but majority of these phyla have different sets of OTUs demonstrating opposite trends (increasing and decreasing during the soil formation). 6. The microbiome of carbonate soils in all cases of benchmark and newly formed soils is familiar because of the uniform zonal environmental conditions in geochemically specific landscape of Izhora upland. It should be taken into account that specialists will elaborate new technologies of soil reclamation in case of humid environment chronologically based on carbonate-enriched heap mining materials.