Microbiomes of different ages in Rendzic Leptosols in the Crimean Peninsula

Study sites and sample collection

All sites are represented by the Rendzic Leptosols located in the first and second mountain ridges in the forest-steppe zone of Crimean Peninsula. The climate of this zone is more humid than in northern part of the Peninsula. The annual precipitation rate is about 380–500 mm per year while the evaporation rate is 750–850 mm. Annual average temperature is +20–22 °C. The depth of soil freezing is no more than 20 cm. Overall, the climate of the investigated area is very close to Mediterranean one. The heights of the relief range between 300–750 m, while topography of the territory is strongly affected by composition and texture of limestones. Limestones are presented by sedimentary rocks strongly affected by karst and denudation processes. Originally, limestone surface was not covered by any other quaternary sediment and this fact provides the possibility for the soil to be formed according to the model of primary soil formation. Thus, all sites are comparable in terms of pedogenesis conditions. Meanwhile, all sites comprise different chronosequence stages, which originated from anthropogenic exploitation of mines for construction and other processes in different historical periods. Age of each chronosequence stage was confirmed by historical documents (Lisetskii & Ergina, 2010). Benchmark site K3 was presented by native brown soil, formed around the Holocene. Site K1 with the oldest anthropogenic impact is located in the 700-year-old territory of medieval fortress city Eski-Kermen, which was destroyed at the end of the 14th century. Near K3 site is site K2, representing 75-year-old WWII trenches in Holmovka village. Site K6 is an overgrown quarry in the north of Belogorsky district with gravel-sandy textured carbonate containing heaps, which was reclaimed approximately 50 years ago. All soil profiles are Leptosols of various thicknesses; the thickness of the humus horizon and the degree of weathering of the fine earth soil increased with age. Samples were collected in the summer of 2018. They were taken from each soil profile for each horizon in 5 replicates. Quantity of horizons differed through sites due to differences in soil profiles: O, AY and C from the K1, AY and C from K2 and K3, AY from K6. The coordinates of the K1 site were 44°36.554 N, 33°44.376 E; K2 and K3 sites 44°39.171 N, 33°44.968 E; K6 site 45°07.644 N, 34°35.537 E (Fig. 1). All soil samples were acquired with the approval of V.I. Vernadsky Crimean Federal University.

Soils for routine analyses were ground and passed through a 2 mm sieve; the large root debris was removed manually. The main agrochemical parameters were measured: P₂O₅ and K₂O by the Machigin method (GOST 26205-91, 1991), pH (GOST 26213-91, 1991) and total nitrogen (GOST 26107-84, 1984). Total organic carbon (TOC) was determined using a CHN analyser Leco CHN-628 (Leco Corporation, USA) in Research Park of St Petersburg State University.

DNA isolation, real-time and 16S rDNA library preparation

For the microbiome analysis, five replicate soil samples were collected from each horizon from each site (40 total samples). From each sample, total DNA was isolated from 0.5 gram of soil using the NucleoSpin®Soil Kit (Macherey-Nagel GmbH & Co. KG, Germany) using a combination of SL1+SX buffers, recommended for soils with low organic content (Lazarevic et al., 2013). Samples were mechanically disintegrated using a Precellys 24 homogeniser (Bertin Technologies, France). The quality of isolation was tested by gel electrophoresis in 1% agarose gel (0.5 × TAE buffer). DNA concentrations were measured at 260nm using SPECTROstar Nano (BMG LABTECH, Ortenberg, Germany). The final DNA concentration was, on average, 50 ng/μL.

Quantitative PCR (qPCR) was conducted for two groups of organisms: bacteria and archaea as previously described in Gladkov et al. (2019). Each sample, including standards, was analysed in triplicate. The mean values with standard errors were calculated for replicates of both PCR and DNA samples. After processing, the results were expressed as the common logarithm of the number of ribosomal operons per 1 g soil.

Construction and sequencing of the 16S rRNA amplicon libraries was conducted using an Illumina MiSeq (Illumina, Inc, USA) at the Centre for Genomic Technologies, Proteomics and Cell Biology (ARRIAM, Russia) as described in Gladkov et al. (2019).

Data processing

Amplicon libraries of the 16S rRNA gene were processed using packages in the R (R Core Team, 2018) and QIIME2 (Bolyen et al., 2019) software environments. RStudio Team (2016) was used as the development environment for R. Trimming, combining sequences into phylotypes and subsequent processing was performed through dada2 package (Nearing et al., 2018), which provides more reproducible and accurate results due to the use of denoising algorithms rather than clustering of phylotypes, in contrast to more classical approaches (Callahan et al., 2016). The taxonomic affiliation of phylotypes was determined using the RDP classifier (Wang et al., 2007) based on Silva 132 (Quast et al., 2013). The phylogenetic tree was built in the QIIME2 software environment using the SEPP package (Janssen et al., 2018). For some analyses, data were normalised by phyloseq (McMurdie & Holmes, 2013) using the rarefaction algorithm according to the sample with the smallest number of readings, and were stabilised by variation through the Deseq2 package (Love, Huber & Anders, 2014) to compare the relative abundances of phylotypes in the samples. For the analysis of alpha diversity, the following indices were used: the observed OTU, Shannon (Shannon & Weaver, 1949), the inverse Simpson (Simpson, 1949) and Faith phylogenetic diversity (Faith, 1992). Significance of mean differences was calculated by the Mann–Whitney test (Mann & Whitney, 1947). For the analysis of beta diversity, communities were compared using the construction of their dissimilarity matrix using the weighted UniFrac, unweighted UniFrac (Lozupone & Knight, 2005) and Bray-Curtis algorithms (Bray & Curtis, 1957). When visualising the data on beta diversity, the dimensions of the dissimilarity matrices were reduced using NMDS (Kruskal, 1964). The significance of sample separation in the analysis of beta diversity was assessed by PERMANOVA (Anderson, 2017) in the form of the adonis2 test as part of the vegan package (Oksanen et al., 2019). To analyse the variation of beta diversity by soil chemical parameters, the constrained correspondence analysis (CCA) was used (Ter Braak, 1986; Palmer, 1993; McCune, 1997). To assess the possible multicollinearity of the CCA model, generalised variance-inflation factors for linear models were used (Fox & Monette, 1992; Fox, 1997). The CCA function and reliability analyses of the model were conducted using the vegan package. To estimate the significance of differences between phylotypes, previously normalised data were processed using the Wald test, with a Benjamin-Hotchberg false discovery rate (FDR) correction in the DEseq2 package (Benjamini & Hochberg, 1995).

The R packages phyloseq, ggpubr (Kassambara, 2019), picante (Kembel et al., 2010), ggforce (Pedersen, 2019), tidyverse (Wickham et al., 2019), ggtree (Yu et al., 2018), ampvis2 (Andersen et al., 2018) were used for post-processing and visualisation of the obtained data.

Soil chemical parameters

All soils demonstrated alkalinity (from 8.2 to 7.6) and a high content of carbonates (4.8–45.6%), which is typical of Rendzic Leptosols. For the K1 and K2 sites, pH and carbonates decreased towards the upper horizons (topsoil), due to leaching processes. Carbonate content in the C horizon at the K3 site (4.8%) was lower than in the AY horizon (28.57%) because most of the carbonates are immobilised in the skeleton of the soil. K1 was the only site with an O horizon in a soil profile; this type of horizon is formed by herbage without grazing. Therefore, it had the highest amount of total organic carbon (TOC) and nitrogen content. Leptosol in site K6 was of a slightly alkaline pH (7.7) and had significant reserves of potassium (1,110 mg/kg) and phosphorus (285 mg/kg), caused by the use of the surrounding territory by the locals of the Vishennoe village to dispose of household waste (Table 1).

Quantitative PCR

Quantitative PCR showed that the number of bacteria ribosomal operons per 1g of soil was high across all sites and horizons (Fig. 2). The archaea operon count varied by horizon, but for K1 and K2 sites, it increased towards the lower horizons.

Initial quality control and phyla composition

After initial processing of 40 amplicon libraries of 16S rRNA genes, three samples were excluded from the subsequent analysis due to their poor agreement with a rarefaction curve (Fig. S1). All data is available at SRA database (SRA Toolkit Development Team, 2020) under BioProject ID PRJNA645404. The final output of the 16S rRNA gene library sequencing included 37 samples with a total of 1145454 reads. The minimum number of detected reads were 13,925, maximum – 41,384 and average read number – 30,958.22. A total of 12,311 OTUs were observed: 11,705 (95%) OTUs were assigned to the kingdom level, 11,026 (89.56%) –phylum level, 10,814 (87.84%) –class level, 9406 (76.4%) –order level, 7800 (63.36%) –family level, 3993 (32.43%) –genus level and 277 (2.25%) –species level.

The most abundant phyla across all samples were Actinobacteria, Proteobacteria, Acidobacteria, Bacteroidetes, Thaumarchaeota, Planctomycetes, Verrucomicrobia, Firmicutes and Chloroflexi (Fig. 3). Site K1 was most clearly distinct from other sites by phyla composition, the most drastic difference being the almost complete absence of Firmicutes representatives. Some phyla demonstrated shifts in abundance correlating with the soil horizon: Bacteroidetes and Proteobacteria were more abundant in topsoil horizons, while Thaumarchaeota, Acidobacteria and Verrucomicrobia were more abundant in lower horizons. This observation is consistent with qPCR data. Bacterial ribosomal count was approximately equal between horizons probably because different bacteria groups shift their abundance in opposing directions through the horizons; and Thaumarchaeota, being the dominant archaeal phyla, was responsible for total archaea increase in lower horizons.

On a family level, the most abundant taxa were Nitrososphaeraceae (Thaumarchaeota), Chitinophagaceae and Microscillaceae (Bacteroidetes), 67-14 and Micromonosporaceae (Actinobacteria), Xanthobacteriaceae and Burkhorderiaceae (Proteobacteria) and Pyrinomonadaceae (Acidobacteria) (Fig. S2). In K1 site samples, phylotypes from Rubrobacteriaceae and Bacillales were less abundant than in other samples and Solirubrobacteriaceae were more abundant. Sphingomonadaceae were more abundant in topsoil. Xiphinematobacteriaceae were more abundant in deeper AC and C horizons.

Alpha diversity

All alpha diversity indices revealed the higher horizons demonstrated a tendency towards higher diversity (Fig. 4). The maximum observed number of OTUs was detected in K1-O and K6-AY, and the minimum in K3-C. In the K1 and K3 sites, observed OTU significantly decreased toward the lower horizon. AY horizons across all sites had comparable OTU numbers. The Faith index, which demonstrates phylogenetic distance (PD), was evenly distributed between samples, with no apparent maximum or minimum. However, it also showed separation of samples by horizon. The Shannon index evaluates diversity, particularly evenness, with respect to minor taxa, while the inverted Simpson index takes into account more abundant taxa. Using the Shannon index, K1-O was similar in diversity to K6-AY but was different according to the inverted Simpson index. In general, the inverted Simpson index shows that K1-O was most diverse, while samples from K2, K3 and K6 sites show significant, but slight, separation from each other. Furthermore, K6-AY was closer in diversity to K1-O and K1-AY than K2 and K3 sites by the Shannon index. In summary, all diversity indices to varying extents show a separation of samples by soil horizon, alongside with secluded position of samples from K1 and K6 sites.

Beta-diversity and CCA

Beta diversity demonstrated two clear trends, coinciding with the axes (Fig. 5). Along the “Y” axis samples tended to line up according to the soil horizons. Along the “X” axis, samples were divided into “site” groups: Bray-Curtis and UniFrac algorithms showed that one group included all samples from the K1 site, the second group included the only sample from K6 site and the third, all samples from both K2 and K3 sites. According to the weighted UniFrac algorithm, samples from the K6 site group together with samples from K2 and K3 sites, which is consistent with results of the inverted Simpson index.

PERMANOVA showed that soil horizon had the maximum coefficient of determination (Table 2). The next factor was the sampling site. All the soil agrochemical parameters, except for carbonates (C_carb), demonstrated similar significance, but with low coefficient of determination values. PERMANOVA nested by horizon showed that all agrochemical factors, including C_carb, became significant (Table S1).

The model of CCA performed for agrochemical factors is statistically significant, although it demonstrated that these factors could not explain the discrepancy between sample sites (Fig. 6). However, they explained soil stratification into horizons. The test on the variance inflation factor showed all agrochemical factors, including pH, demonstrated multicollinearity. A combination of CCA and PERMANOVA confirm that variability between soil horizons was associated with agrochemical factors.

K1/K3 phylotype comparisons

Previous analyses concluded that microbiomes across all sites separate by soil horizon, but also that microbiomes from the K1 site are more distinct from other sites. To assess more precise differences in microbiome composition between sites, we visualised significant shifts of phylotype abundance in AY and AC/C horizons between the K1 and K3 sites (Fig. 7). Despite the major trend of microbiome differences between soil horizons, our analysis shows that the reactive component of the soil microbiome shifted together in both soil horizons between different soil sites. Firmicutes, in particular Planococcaceae and B. longiquaestium, increased in K3; Actinobacteria (Solirubrobacter, Gaiella, 67-14, Microlunatus, Ilumatobacteraceae) mostly increased in K1, except for Rubrobacter; Proteobacteria (Deltaproteobacteria, Bradyrhizobium, Xanthobacteriaceae, Rhodoplanes, Pedomicrobium, Reyranella, Geminicoccaceae, Burkhordeliaceae, MND1, Steroidobacter uvarum) were more abundant in K1. Representatives of Verrucomicrobia (Xiphinematobacter, Udaeobacter), Thaumarchaeota (Nitrososphaeraceae) and Acidobacteria (NA, RB41) varied in both K1 and K3 sites. Variation of Thaumarchaeota in both K1 and K3, the growth of which depends on nitrogen content, supports the earlier conclusions that nitrogen content doesn’t explain site differences. However, the K1 site was abundant in Actinobacteria and Proteobacteria related phylotypes.

K2/K3 phylotype comparisons

Microbiomes of the AY and C horizons from two sites in Holmovka village (K2 and K3) were the closest to each other on the beta diversity plots. These data are supported by log2FoldChange values for the 30 of the most abundant phylotypes of both horizons between sites (Table S2). Almost half of these phylotype changes were not significant. The greatest differences in topsoil (more than 10 times more at K3-AY than K2-AY) were for Seq13 (Oxyphotobacteria from Cyanobacteria), Seq101 and Seq136 (Planococcaceae from Firmicutes), Seq322 (Chitinophagaceae from Bacteroidetes) and Seq339 (Romboutsia from Firmicutes). For the deeper horizon, the only phylotype matching these conditions was Seq445 (Adhaeribacter from Bacteroidetes).

K6/K3/K1 phylotype comparisons

To assess the specificity of microbiome composition of the Leptosol at the K6 site, similar to K2 and K3 sites, we estimated shifts in abundance by calculating log2FoldChange values for 30 phylotypes for K6-AY/K3-AY and K6-AY/K1-AY pairs (Table S3). All log2FoldChange values were significant, except for the only phylotype in the K6-AY/K3-AY pair. Eleven phylotypes appeared in both pairs of comparisons, and most of them were more abundant in K6-AY: Seq20 and Seq161 (Nitrososphaeraceae from Thaumarchaeota), Seq11 and Seq119 (RB41 from Acidobacteria), Seq94 (Subgroup 6 from Acidobacteria), Seq53 (Chitinophagaceae from Bacteroidetes) and Seq165 (Oxyphotobacteria from Cyanobacteria). However, a lot of other phylotypes were underrepresented in K6 compared to the other two sites. In comparison with K3-AY, the K6-AY site contained more than 10 times fewer of the following phylotypes: Seq13 (Oxyphotobacteria from Cyanobacteria), Seq5 (Candidatus_Xiphinematobacter from Verrucomicrobia), Seq37 (Chitinophagaceae from Bacteroidetes), Seq101 and Seq136 (Planococcaceae from Firmicutes), Seq60 (Aridibacter famidurans from Acidobacteria) and Seq34 (Thermoleophilia from Actinobacteria). Compared to K1-AY, the K6-AY site contained more than 10 times less of the following phylotypes: Seq6 and Seq36 (Thermoleophilia from Actinobacteria), Seq33 (Nitrososphaeraceae from Thaumarchaeota), Seq25 (Microlunatus from Actinobacteria) and Seq49 (MND1 from Proteobacteria). K6-AY was more abundant than K1-AY by Seq3 (Bacillus longiquaesitum from Firmicutes) and Seq32 (Candidatus _Nitrososphaera from Thaumarchaeota). These differences show that topsoil microbiomes of all sites were composed of similar major phylotypes, including both oligo- and copiotrophic taxa, which shifted between sites regardless of their trophic group. These data are consistent with the observation that the changing soil chemical parameters have not explained the beta diversity observed between sites.