Diversity of Arbuscular Mycorrhizal Fungi in Distinct Ecosystems of the North Caucasus, a Temperate Biodiversity Hotspot

1. Introduction

One of the key factors of species diversity conservation is the existence of biodiversity hotspots – regions with extremely high levels of species richness, endemism, and number of threatened species [1,2]. The conception of hotspots was suggested originally for tropical forests [3], but with intensive investigation of global biota many various hotspots have been discovered on the Earth. There are 36 biodiversity hotspots according to the last revision [4]. The hotspots conception has been subjected to considerable criticism (summarized in [2]), mostly because the majority of studies in the field of hotspots biodiversity have only focused on plant and vertebrate species, whereas the diversity of invertebrates and fungi has not been studied enough [5]. The modern approach to biodiversity conservation should be more complex and includedifferent aspects, like interactions between various organisms [2]. Arbuscular mycorrhizal fungi (AMF) can play an important role in addressing biodiversity conservation. AMF are microscopic soil fungi that form the most common type of plant-microbial symbiosis. A total of 71 % of terrestrial plants are characterized by obligate symbiotic relationships with AMF and only 7 % of plants form inconsistent “nonmycorrhizal – mycorrhizal” associations [6]. Plants receive a number of benefits from this interaction: fungi improve the uptake of phosphorus, other nutrients and water, increase resistance to pathogens, and protect plants from toxic heavy metals [7]; Edlinger et al., 2022; Kakouridis et al., 2022.

Research of AMF diversity in some tropical biodiversity hotspots has revealed a significant richness of species in these areas. For instance, 49 AMF species were identified in Cerrado hotspot, Brazil [8]; 58 AMF species were found in the rhizosphere of plants from the Atlantic Forest [9]; 72 species were discovered on a tropical mountain in Brazil [10]. While some research has been carried out on tropical AMF, few studies have investigated AMF diversity in temperate hotspots.

One of the richest temperate regions on the Earth is the Caucasus [11]. Since 2000 it has been considered as a biodiversity hotspot [1]. The Caucasus is a region on the isthmus between the Black and the Caspian seas, on the border between Europe and Asia. The Caucasus ecoregion is shared by Russia (Karachay-Cherkess Republic, Republic of Dagestan, Republic of Ingushetia, Kabardino-Balkarian Republic, Republic of North Ossetia–Alania, Chechen Republic), Georgia, Armenia, Azerbaijan, Turkey, and Iran. The Karachay-Cherkess Republic located on the northern macro slope of the Caucasus Mountains [12] is a promising region for the study of biodiversity. It is characterized by a complex relief, a variety of natural obstacles, frequent changes of ecotopes, mosaic vegetation cover, and phytocenoses diversity in the mountainous southern part of the region. The main waterway is the Kuban river, which has only left tributaries descending from the Caucasus Mountains (such as the Teberda river). Investigations on AMF diversity in biotopes of river valleys, forests, alpine and subalpine meadows of the Karachay-Cherkessia are very limited. While previous studies of fungal biodiversity in Caucasus have been carried out at the family [13] or the genus level [14] by molecular genetic methods, no studies have addressed AMF biodiversity in Caucasus at the species level. Twelve AMF species were found in soil-root samples from southern Georgia [15]. But this study employed only morphology-based identification, whereas the molecular-based methods allows to identify more taxa from the same soil and roots samples than the morphological approach [16,17].

To assess molecular diversity of AMF, researchers use the DNA-metabarcoding approach, including high-throughput sequencing of the DNA barcode from soil and root samples, followed by bioinformatics analysis. The critical point for fungal metabarcoding is the length of obtained sequences, therefore high-throughput sequencing platforms able to produce long reads are used. Before 2020 mainly pyrosequencing (454-Roche) was employed [16,18,19,20], now Illumina MiSeq prevails in fungal metabarcoding studies [17,21,22,23]. Another promising platform is PacBio, which can obtain a sequence of 2.5 kb of rDNA in one read [24].

There are no consensus metabarcoding markers for AMF, but almost all currently employed barcodes are parts of the rDNA cistron. In the last decade internal transcribed spacers (ITSs) of rDNA have become the most common universal fungal metabarcoding marker [25]. Fragments of SSU (small sub-units, or 18S) rRNA gene [16,20,22]; ITS [19,26]; and fragments of LSU (large subunit, or 28S) rRNA gene [18,27] are frequently used as barcodes for AMF. For amplification of SSU rRNA gene fragments, AMF-specific primers pairs are being developed (for example, [28]). In case the goal of the study is to observe not only AMF, but the diversity of various fungi classes, ITS and D1/D2 regions are used [18,19]. For identification based on ITSs, universal eucaryotic primers or a combination of fungi-specific primers with universal primers are used [19].

Besides rDNA other DNA markers are seldom used as barcodes. The pilot experiment showed that glomalin, a homolog of a heat shock protein, is a good marker for AMF identification [29]. However, the major drawback of introduction of new barcodes is a lack of reference data. For successful identification of species, a reference database should contain a representative number of correctly identified taxa. Almost 350 species of AMF have been correctly described according to the Botanical Nomenclature Code [30]. The biggest public database GenBank NCBI [31] contains glomalin data only for 24 species. In addition to glomalin, the following genes and regions were used as barcodes: 1) the mitochondrial cytochrome c oxidase 1, COX1 or COI [32,33]; 2) LSU of mitochondrial rRNA [34]; 3) I and II subunits of RNA polymerase II [35]; 4) H+-ATPase [36]; 5) β-tubulin [37]; 6) actin; 7) elongation factor 1-alpha, EF 1-α; 8) phosphate transporters; 9) RNA polymerase II subunits (RPB1 and less RPB2) [38,39]. Most of these markers have not been widely used due to the low nucleotide polymorphism between closely related AMF species. Perhaps the rpb1 gene can become the second barcode in quality after 18S-ITS1-5.8S-ITS2-28S region, as rpb1 was successfully used for identification of AMF species [39,40].

The systematics of AMF is quite complex and controversial. There are many synonyms; many genera were taxonomically revised, but after the revisions not all AMF species names were effectively published according to the Botanical Nomenclature Code. Many taxa found in soil samples are described only by the DNA barcode and have no proper scientific names. AMF form the Glomeromycota phylum and are classified into one class, Glomeromycetes, and four orders [Schüßler, 2023]. But, according to Castillo et al. [42], AMF form the Glomeromycota phylum and are classified into three classes and five orders. According to the database curated by SYMPLANTA (Symbiosis and Plant-Microbe Association Research Laboratory), there are four orders, 12 families, 44genera and 443 species of AMF fungi [30]. Various genera are represented in the GenBank NCBI unevenly. The database contains the majority of species from Claroideoglomus, Dominikia, Septoglomus, Corymbiglomus, Diversispora genera. Poorly represented genera include Glomus sensu lato, Scutellospora, Paraglomus, and Pacispora. Thus, the GenBank database is not complete and lacks many species. Due to the problem of precise identification of AMF species, dozens of virtual taxa were obtained [17,20,21,43]. The alternatives for GenBank NCBI database are the specialized reference databases. One of the main curated databases is MaarjAM that includes mainly SSU rRNA gene-based barcodes and checked metadata for each entry [44]. This database is regularly updated and now contains data for different barcodes [45]. We can assume that the most accurate assessment of biodiversity of AMF (Glomeromycetes) and other fungal classes is their identification using rDNA ITSs as universal barcodes.

One of the most important directions of modern studies is to identify the key factors causingAMF biodiversity. Among them are: temperature conditions, agrochemical properties of soils (including the pH level), usage of meadows as grassland, changes of plant communities in the process of succession [46–49]. According to Horn et al. [50], biotic factors have a more significant impact on the composition of AMF communities in comparison with abiotic factors. The development of succession from a meadow to a forest accompanied by an increase in the proportion of perennial plants in the ecosystem may be the reason for a decrease in the abundance of AMF [51,52]. According to Torrecillas et al., annual plants had a higher AMF diversity [53]. However, no convincing unambiguous evidence has been obtained for the presence of a key factor determining the high species diversity of AMF. Another factor that is consideredas an important one is an altitude. For example, a subalpine meadow (>2,500 m above sea level) in the Head Mountain, Inner Mongolia Autonomous Region, was shown to have the greatest diversity of AMF [49]. On the other hand, in the Southern-Central zone of Chile, the “natural grassland” ecosystem had the greatest diversity of AMF species at a significantly lower altitude (<500 m above sea level) [42]. Perhaps, there are no direct correlations between altitude and AMF biodiversity [50,54–57]. Further elucidation of the factors contributing to the increase in AMF species diversity is required. Thus, it is essential to test the supposition that the biodiversity of AM fungi may be primarily related to such factors as soil pH, phosphorus content, altitude above sea level, and plant species richness. This study aim to detect the AMF hotspots and define the AMF biodiversity in rhizospheric soils of different biotopes in the North Caucasus (such as a subalpine meadow, a forest and a river valley). Research objectives are: 1) to sequence fungal DNA samples from the rhizospheric soil using Illumina MiSeq and ITS1 and ITS2 regions as barcodes; 2) to assess the species diversity of AM fungi; 3) to assess the key relationships between parameters of the biodiversity of AM fungi and plant species richness, altitude for trial plots, and agrochemical soil parameters.

2. Materials and Methods

2.1. The Characteristics of the Stationary Trial Plots

The characteristics of Stationary Trial Plots (STPs, study areas, sample plot), as well as the soils in which the species diversity of AMF was determined, are presented in Table 1 and Table S1. The following STPs were analyzed in different altitudes: 1) in subalpine meadow biotopes – STP 1 (subalpine meadow-4, Malaya Hatipara ridge, 2,437 m above sea level), STP 3 (subalpine meadow-3, Malaya Hatipara ridge, 2,401 m above sea level), and STP 4 (subalpine meadow-2, Malaya Hatipara ridge, 2,186 m above sea level); 2) in forest biotopes – STP 7 (fir forest-3, Malaya Hatipara mountain, 1,900 m above sea level), STP 8 (pine forest-3, Malaya Hatipara mountain, 1,890 m above sea level), and STP 9 (mixed forest near the Bolshaya Hatipara river, Bolshaya Hatipara mountain, 1,507 m above sea level); 3) in river valley biotopes – STP 11 (grassland in the valley of the Teberda river, Teberda town, 1,342 m above sea level), STP 12 (grassland in the valley of the Teberda river, the border of New Teberda village, 1,026 m above sea level), and STP 13 (grassland in the valley of the Kuban river, Ordzhonikidzevsky village, 795 m above sea level). The trial area was 10 m x 10 m. A botanical survey of test sites, STPs located in undisturbed ecosystems of the North Caucasus, the Teberdinsky National Park and adjacent territories, was carried out using A.S. Zernov’s determinants ([12,58]; see Table S2). In each STP, the soil samples were collected for molecular genetic analysis of fungal diversity from the upper 0-5 cm horizon of the soil without litter according to [14], and the samples of the upper horizon of the soil without litter were collected for agrochemical analysis according to [59].The total phosphorus content (P_total, %) by the Ginzburg&Shcheglova's method, the content of phosphorus available for plant nutrition (Pi, mg/kg) by the Truog method, the total nitrogen content (N_total, %) by the Kjeldahl method, pH_KCl, the sum of fractions <0.01mm, soil type- by common methods were determined ([60]; see Table S1).

2.2. Molecular Genetic Identification of Fungi

2.2.1. Sampling and Molecular Analysis

For the molecular genetic identification of AMF, a sampling technique for Illumina MiSeq sequencing optimized by the authors was used. In total, 10 samples were taken from each STP. The samples were frozen and transported in 50 ml tubes in liquid nitrogen. Samples of rhizospheric soil containing mycelium and AMF spores were taken for identification (verification by microscopy). To isolate DNA, 0.2 g of frozen soil and 0.5 g of granite sand were taken into a 2 ml test Eppendorf tube for mechanical grinding (degradation) of the material. After that, 700 µl of СТАB buffer was added to the heated test tube with soil, (2% СТАB; 1,4M NaCl; 20mM EDTA; 100mM Tris-HCl pH=8.0). Initially, the material was homogenized in a buffer with granite sand for 15 minutes on a vortex. After that, the test tubes were shaken every 15 minutes by laboratory vortex shaker (Biosan, Latvija) for 1 min during 1 hour at +65 °C. After thermal, chemical and mechanical treatment, the samples were centrifuged for 5 minutes and the supernatant was transferred to new test tubes. DNA was additionally washed off soil residues by 500 mkl of water by shaking soil with water in the vortex for 5 minutes. After additional centrifugation, the second supernatant was combined with the first one. Purification of the obtained DNA was carried out by double extraction with an equal volume of chloroform. Each time after centrifugation (10 min at 14000 rpm, Eppendorf, Germany), a DNA supernatant was selected. Then the DNA was deposited with 2/3 V isopropanol with 0.4 M NaCl, and the precipitate was washed with 70% ethyl alcohol, the precipitate was dried (until the alcohol was completely removed) for 3 minutes and then dissolved in water [61,62]. After that DNA was purified with the Qiagen Gel Extraction Kit (Germany). Purification on columns or with silicon oxide after isolation of a DNA fragment on an agarose gel after electrophoresis does not allow to completely get rid of various PCR inhibiting impurities. Therefore, an additional measure could be the dilution of DNA samples before amplification by 100 times with water (accompanied by diligent mixing) in order to reduce the effect on PCR of inhibitors from the soil, including humic acids.

Purified DNA was used for nested PCR marker regions ITS1 and ITS2 with universal primers. The primers were synthesized in Evrogen, Russia with adapter 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3’ for direct primers and adapter 5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3’ for reverse primers for Illumina MiSeq, such as: primer ITS5 (5’-GGAAGTAAAAGTCGTAACAAGG-3’) and reverse primer ITS-2RK (5’-CGTTCAAAGATTCGATGATTCAC-3’) modified by the authors for amplification of the ITS1 region, as well as primers ITS3 (5’-GCATCGATGAAGAACGCAGC-3’) and reverse ITS4 (5’-TCCTCCGCTTATTGATATGC-3’) for amplification of the ITS2 region. PCR: 1) PCR for the long fragment (ITS1 + ITS2 regions) with primers ITS5 and ITS4; thermal cycling conditions were: 5 min initial denaturation and polymerase activation at 95 °C; 35 cycles of 20 s denaturation at 95 °C, 20 s annealing at 58 °C and 40 s elongation at 72 °C; and final elongation at 72 °C for 10 min; 2)verification of PCR with gel electrophoresis, DNA dilution by 1000 times; 3) 35 PCR cycles into short fragments (separately for ITS1 and ITS2 regions); repetition of the program with primers for nested PCR; 4) quality control with gel electrophoresis, PCR products for ITS1 and ITS2 were combined for each sample and purified on magnetic particles AMРure XP (Beckman Coulter, USA).

Amplicon libraries were sequenced on an Illumina MiSeq device using a set of reagents MiSeq® Reagent Kit v3 (600-cycle) with two-way reading (2x300 bp) ("Illumina, Inc.", USA). The identified sequences were processed using Illumina software ("Illumina, Inc."). As a result of sequencing on the Illumina MiSeq platform, FASTQ sequences were obtained from forward and reverse primers. This format includes sequence records and reading quality indicators (quality score) for each nucleotide position.

2.2.2. Bioinformatics

We uploaded the Illumina data (FASTQ sequences) to the NCBI database; it is now available in the base as the PRJNA646244 bioproject [Bioproject PRJNA646244, 2023]. Data processing was performed with USEARCH [63]. Forward and reverse reads were merged using the command “fastq_mergepairs” with sequences shorter than 120 bp removed. The resulting reads were pooled into a single FASTQ file. Further processing to OTUs with a 97% radius was carried out in two steps. First, sequences from the MaarjAM database [44] were used as OTU centroids, to which the obtained sequences in FASTQ format were mapped. Secondly, the sequences that have not been mapped to the MaarjAM database (with 97% threshold) were processed to OTUs using the UPARSE algorithm [64], where the most abundant sequences are selected as centroids. During this step regions corresponding to primers were removed, and sequences were quality filtered with the value of “fastq_maxee” parameter was equal to 1. Chimeric sequences were removed during OTU construction using the “cluster_otus” command, and residual singletons were removed. Both types of centroids were used to construct the OTU table. The high genetic polymorphism observed in the internal transcribed spacers (ITSs) of AMF [Thiéry et al., 2012] makes it difficult to identify many AMF species accurately [GenBank NCBI, 2023]. Consequently, an increasing number of virtual taxa (VT) are being added to databases [Peña-Venegas et al., 2019; Lara-Pérez et al., 2020; Rasmussen et al., 2022]. For effective assessment of AMF diversity, Operational Taxonomic Units (OTUs) is recommended instead of Amplicon Sequence Variants (ASVs) or Zero-radius OTUs (ZOTUs) for the ITS1 and ITS2 AMF regions.

For taxonomic annotation of the newly obtained OTUs we used the MaarjAM database, as well as the established local reference database consisting of sequences of fungi of the Glomeromycetes class selected from the NCBI database. The annotation was performed using the “usearch_global” command. The coding regions of rRNA were trimmed, leaving the spacer sequences for further processing. Primary clustering was performed by constructing an alignment-free tree using the Neighbor Joining (NJ) method in the paHMM program [65]. Based on this tree, sequence clusters that could be properly aligned were selected. These alignments were supplemented with reference sequences annotated at the species level and were used to construct trees using the Maximum Likelihood (ML) method in the IQ-TREE [66]. Fungi annotation at the phylum level was based on the UTAX reference dataset version 27.10.2022 [67].

The general algorithm for further identification is shown in Figure S1. The alignment of ITS sequences by genera was most effective, since the inter-genus variability in ITS for arbuscular mycorrhiza (AM) fungi is too high and most often did not allow for reliable alignment of sequences, and this is the most significant factor for the construction of reliable phylogenetic trees. With the help of the Mega 7 software package [68], sequences combined by genera into separate files were aligned automatically, and then a manual revision of the alignment was carried out. After manual refinement of the alignment and the BLAST procedure in the NCBI GenBank, some sequences were rejected as not belonging to this genus (rarely, to a taxon of a higher rank), due to the inability to align them within the ITS regions under study. During manual alignment, species-specific mutations were also identified, the analysis of which provided additional information about the similarity of sequences. In some cases, this was the only way to separate sequences of two close species, which is most important in particular for the identification of species from the Rhizophagus genus. After that, phylogenetic trees were reconstructed in the Mega 7 software package. The resulting trees were used to refine the identification at the species level. Sequences clustered together to form a separate clade were identified as sequences belonging to the same taxon. If the resulting clade was found to include reference sequences from the GenBank belonging to the same species, such sequences were assigned to a particular species. To clarify the species affiliation we also used information on the length of phylogenetic branches as well as p-distances calculated by us using the NCBI database for reference sequences (Table S3), bootstrap indexes and the presence or absence of species-specific nucleotide substitutions (e.g. [69]). In complex cases, individual sequences were also identified taking into account the closest sequences based on the results of manual BLAST analysis in the NCBI GenBank, paying attention primarily to the year of obtaining the sequence presented in the GenBank. If the obtained sequences were grouped into a clade, but not grouped with reference sequences, then they were assigned to a virtual taxon (VT) of the species level. If VT was not clustered on a tree with other sequences, forming a well-separated clade and had low p-distance compared to others, then it was identified as VT of the genus level or higher. If such VT consisted of multiple sequences, they were excluded from the analysis to avoid biases in biodiversity estimation, as it was hard to determine the exact number of species they represented.

2.2.3. Assessment of Biodiversity Indices

To assess the taxonomic diversity, Hill numbers were calculated using the iNEXT package [70] for R environment [71]. Species Richness (diversity order q=0), Shannon diversity (q=1) and Simpson diversity (q=2) were calculated. To overcome the limitations associated with different sample sizes, Hill numbers were also calculated, extrapolated to a sample tending to infinity (estimated asymptotic diversities). A 95% confidence interval was calculated for each extrapolated index.

2.2.4. Statistics

Statistical analysis was processed in the R language environment 4.1.0 [71]. The data, OTUs for fungal classes, were normalized by the total-sum scaling method. The analysis of variance (ANOVA) was used to assess the statistical significance (P<0.05) of the differences. Principal Component Analyses (PCA) was applied using the pcaMethods package [72]. For unsupervised dimension reduction multidimensional scaling (MDS) was applied. Spearman’s distances (1-ρ) and MDS were used with stats [71]. Analysis of similarities (ANOSIM) [Warton et al., 2012] was made with the vegan package [Oksanen et al., 2022], and Spearman's distance (1-rho) was used. Hierarchical clustering was carried out using the dendextend package [73]. Spearman’s distances and Ward method for cluster agglomeration were used. The packages of ggplot2 [74] and Venn Diagram [75] were used to plot the graphs.

3. Results

3.1. Analysis of Fungal Sequences

The total number of fungi sequences in ITS1 and ITS2 regions was >10 million (Table 2). The average depth of reads for STPs of meadows, forests and river valleys (pastures) exceeded 1.1 million, and after filtering (merged reads after length trim) – 376,000 reads (Table 2). The number of AM fungi reads according to ITS1 varied widely for, the average for 9 STPs was 1,774 AMF reads per STP. The results of the Illumina MiSeq sequencing analysis showed that of all the sequences, >62% of the reeds (or >90% of OTUs) belong to fungal sequences. Meanwhile, the number of AMF reads varied more widely for STP in ITS2 compared to ITS1, and the average number of reads for 9 STP was 722 AMF reads per STP in the ITS2 region. Thus, the number of AMF rows in the ITS2 region was on average 2.5 times lower than in the ITS1 region. The highest number of AMF reads was also detected in samples STP 11 and STP 13 related to the river valley, and the lowest in samples STP 7 and STP 8 related to forest biotopes.

It is known that most studies are carried out with the removal of singletons due to the fact that the analyzed objects often make up about 50% of all reads. But the rows of AM fungi represent usually only 1-2% of all reads, so the analysis of singletons can contain important information about their diversity. On one hand, singletons often contain chimeric sequences and can be eliminated at the stage of manual alignment. On the other hand, the exclusion of singletons would lead to a significant loss of true fungal species. The application of public database sequences as OTU centroids has provided a number of advantages. First of all, the use of the same centroids for OTU clustering in different experiments made it possible to compare results in a better way. Secondly, the singletons present in the data were successfully mapped with known centroids, and were not lost at the stages of data processing. Therefore, the use of the MaarjAM database to build OTU was justified, and this approach made it possible to save about 15% of reliable singletons in the final data.

The results showed that for the ITS1 region, the number of OTUs of AMF for the river valley biotopes was 2.5 times higher than for the mountain forest biotopes, and 1.7 times higher than for the subalpine meadow biotopes. The number of OTUs for meadow biotopes was also higher than for forest biotopes (by 1.5 times) according to the ITS1 region. According to the ITS2 region, the number of OTUs for the river valley biotopes was 3.7 times higher than for the forest biotopes, and 1.7 times higher than for the subalpine meadow biotopes. The number of OTUs for meadow biotopes was also higher than for forest biotopes (2.2 times) according to the ITS2 region. Summing up the results of the analysis, we can assume that the biotopes of pastures in the river valley may have the greatest diversity of taxa, while the forest biotopes are characterised by the least diversity. Samples representation in low dimensional space revealed from MDS (Figure S2) and score plots from PCA of OTUs profiles (Figure S3) confirmed the differences between the three studied biotopes.

3.2. Taxonomic Composition of Fungal Phyla and Classes in Soil Samples from the River Valley, Subalpine Meadow and Forest

An assessment of the number and ratio of fungal OTUs at the phylum level was carried out to identify the objective ratios of AMF reads and the general structure of fungal communities in the studied ecotopes (Figure S4). The most common fungal phyla were Ascomycota, Basidiomycota, Glomeromycota, Rozellomycota, Mortierellomycota, Chytridiomycota, Mucoromycota (Figure 1). Kickxellomycota, Olpidiomycota, Zoopagomycota accounted for less that 1% (“Others” at the Figure 1). AM fungi are included in the phylum Glomeromycota, a monophyletic group with one class, Glomeromycetes. Glomeromycota had a significantly higher ratio of OTUs in soil from disturbed ecosystems of pastures in river valley than in soil from forest and subalpine meadow and had a significantly lower ratio of OTUs in soil from forest (P<0.05; Student's t-test). The abundance of ectomycorrhizal fungi (Ascomycota, Basidiomycota) in the analyzed biotopes had no significant differences, except for a reduced level (P<0.05) of Basidiomycota for the river valley. Other fungal phyla also did not simultaneously have significant differences in biotopes, like Glomeromycota.

The following classes of fungi (see Figure 2) were identified in soil samples by ITS1 and ITS2 regions according to the proportion of OTUs for the biotopes of the subalpine meadow, forest biotopes and pasture biotopes in the river valley. A total of 9 fungal classes accounted for at least 2% were found in soil samples from subalpine meadow, forest and river valley (Figure 2). The results show that the same trend can be traced as in the analysis of phyla. With regard to Ascomycetes, we found that Sordariomycetes, Dothideomycetes and Pezizomycetes were more common in the river valley, whereas Leotiomycetes was more common in the forest biotope. AM fungi (Glomeromycetes) were also well represented in soil samples from the river valley, but difference between valley and other biotopes was astonishing.

Thus, the results showed a significant proportion of AMF (Glomeromycetes class) in the total number of sequenced samples attributed to the Fungi kingdom, especially for the river valley. In different STPs there were 0.1-6.9% of all reads and 5-11% of OTUs identified as AMF by ITS1 sequencing. 0.01-1.0% of all reads and 1-5% of OTUs were identified as AMF by ITS2 sequencing. The number of OTUs from Glomeromycetes in the river valley were significantly (P<0.05) larger than in subalpine meadow and OTUs number in subalpine meadow were significantly (P<0.05) more than in the forest (Table 2). Thus the greatest diversity of AMF species can be expected in the biotopes of the river valley pastures (STP 11, STP 12, STP 13).

3.3. Identification of AMF Species in the River Valley, Subalpine Meadow and Forest Biotopes

Intraspecific polymorphism is difficult to assess due to the lack and heterogeneity of data in the NCBI GenBank database. In order to solve this problem we focused on the intra-genus p-distance (Table S3) calculated from the NCBI database when complex (remote from the main species clusters) sequences in the samples. In accordance with the developed algorithm (Figure S1), phylogenetic trees were analyzed for more accurate identification of AMF species, known species were identified, as well as virtual taxa of the species level: an example of Maximum-Likelihood (ML) phylogenetic trees for the genera Ambispora and Acaulospora is shown on Figures S2–S5. The virtual taxa at the species level were defined. Ones made up a significant part of the detected sequences. The construction of phylogenetic trees for all AMF genera, based on the analysis of ITS1 and ITS2 individual regions, resulted in identification of a number of AMF available in the analyzed soil samples. Figures 3 and 4 show heat maps with the frequency of occurrence of identified AMF species by the number of reads for the ITS1 and ITS2 regions.

The analysis of phylogenetic trees revealed that a significant number (>42%) of the detected and identified AMF reads to the species level are virtual taxa. Almost the same proportion (~38%) was made up of AMF reads identified to the level of genus, family and higher rank of taxa (it was not included in the performed analysis because cannot characterize species diversity). The assumption is that at the moment the AMF discovered in the North Caucasus are not fully annotated by genetic methods. Only the construction of trees separately by genera (or by groups of close genera) is supposed to identify the taxa at the species level with a high variability of genetic markers. By constructing phylogenetic trees, the following AMF taxa defined up to the species level were reliably identified in STPs in three different biotopes. The following known species have been identified for ITS1 (see Table S4):

-

subalpine meadow (Acaulospora alpina, Aс. nivalis, Aс. nivalis, Aс. paulinae, Aс. punctata, Aс. viridis, Ambispora gerdemannii, Am. leptoticha, Archaeospora europaea, Ar. trappei, Claroideoglomus claroideum, Cl. walkeri, Diversispora insculpta, Dominikia bernensis, Do. disticha, Glomus indicum, G. macrocarpum, G. tetrastratosum, Otospora bareae, Paraglomus brasilianum, Paraglomus laccatum, Rhizophagus intraradices, Rhizophagus irregularis, Septoglomus constrictum, S. nigrum), or 24 species from 12 genera, as well as 19 virtual taxa at the species level;

-

forest (Acaulospora nivalis, Ac. paulinae, Ac. punctata, Ambispora fennica, Am. gerdemannii, Am. leptoticha, Claroideoglomus claroideum, Diversispora insculpta, Di. slowinskiensis, Di. sporocarpia, Di. spurca, Dominikia bernensis, Glomus indicum, G. macrocarpum, G. tetrastratosum, Otospora bareae, Paraglomus laccatum, Rhizoglomus invermaium, Rhizophagus intraradices, Rhizophagus irregularis, Septoglomus constrictum, S. nigrum), or 22 species from 11 genera, as well as 12 virtual taxa at the species level;

-

river valley (Acaulospora delicata, Ac. paulinae, Ambispora fennica, Am. gerdemannii, Archaeospora europaea, Ar. trappei, Claroideoglomus claroideum, Cl. lamellosum, Cl. walkeri, Diversispora celata, Di. varaderana, Dominikia achra, Do. bernensis, Do. disticha, Funneliformis mosseae, Glomus indicum, G. macrocarpum, G. tetrastratosum, Halonatospora pansihalos, Otospora bareae, Palaeospora spainii, Paraglomus laccatum, Rhizoglomus invermaium, Rhizophagus intraradices, Rhizophagus irregularis, Septoglomus constrictum, S. nigrum, S. viscosum), or 28 species from 15 genera, as well as 24 virtual taxa at the species level.

The list of genera with the number of OTUs identified in the ITS1 region is presented in Table S5. The most common genera (54% OTUs) are highlighted in green; the group of minor genera is highlighted in blue (5% OTUs), which include Ambispora, Otospora, Scutellospora, Palaeospora, Halonatospora.

The following known species have been identified for ITS2 (Table S6):

-

subalpine meadow (Acaulospora alpina, Ac. brasiliensis, Ac. paulinae, Ambispora gerdemannii, Am. leptoticha, Archaeospora trappei, Claroideoglomus claroideum, C. lamellosum, Diversispora varaderana, Dominikia bernensis, Entrophospora infrequens, Glomus bareae, G. indicum, G. macrocarpum, Paraglomus laccatum, Rhizoglomus melanus, Rhizophagus intraradices, Rhizophagus invermaius, Rhizophagus irregularis), or 19 species from 11 genera, as well as 15 virtual taxa at the species level;

-

forest (Acaulospora paulinae, Ambispora gerdemannii, Am. leptoticha, Claroideoglomus claroideum, C. lamellosum, Diversispora slowinskensis, Entrophospora infrequens, Glomus hoi, G. indicum, Paraglomus laccatum, Rhizophagus intraradices, Rhizophagus irregularis, Scutellospora alterata), or 13 species from 9 genera, as well as 7 virtual taxa at the species level;

-

river valley (Acaulospora paulinae, Ac. punctata, Ambispora fennica, Am. gerdemannii, Archaeospora spainiae, Ar. trappei, Cetraspora gilmorei, Claroideoglomus claroideum, Cl. hanlinii, Cl. lamellosum, Dominikia bernensis, Do. difficilevidera, Entrophospora infrequens, Funneliformis mosseae, Glomus hoi, G. indicum, G. macrocarpum, Halonatospora pansihalos, Paraglomus laccatum, Rhizophagus aggregatus, Rhizophagus intraradices, Rhizophagus invermaius, Rhizophagus irregularis, Scutellospora pellucida), or 24 species from 13 genera, as well as 20 virtual taxa at the species level.

The list of genera with the number of OTUs detected for region ITS2 is presented in Table S7. The most common genera (54% OTUs) are highlighted in green; the group of minor genera (5% OTUs) is highlighted in blue. As shown in Tables S5 and S7, ~50% of all AM fungi OTUs are the OTUs belonging to only three genera: Rhizophagus, Dominikia and Glomus, both in regions ITS1 and ITS2. The fungi from these genera were present in OTUs of meadows and river valley, but were absent in some cases in STPs belonging to forest biotopes. Minor genera (5% OTUs) for region ITS2 included Scutellospora, Microkamienskia, Rhizoglomus, Diversispora, Cetraspora, Halonatospora.

The construction of phylogenetic trees indicated that employment of ITS1 guarantied identification of more taxa at the species level. The reason is the specificity of primer annealing on the analyzed genetic markers. The sequence analysis revealed the presence of 560 OTUs identfied up to the species level by regions ITS1 and ITS2 for 9 STPs of three biotopes of the North Caucasus. And 161 OTUs were excluded from the OTUs identified up to the species level due to the lack of clustering in phylogenetic trees constructed for individual AMF genera. 285 OTUs were identified up to 38 species from 16 genera and 27 virtual taxa at the species level by ITS1. 273 OTUs were identified up to 32 species from 15 genera and 25 virtual taxa at the species level by ITS2. Therefore the average number of detected OTUs AMF in the sequences analysis of the ITS1 and ITS2 region was equal (slightly less for ITS2), while the total number of detected reads for the fungi of the Glomeromycetes class was 2,5 times higher by ITS1 than by ITS2 (Table 2).

According to the data of all analyzed STPs the following species were the most common at the culculation of the number of reads for ITS1+ITS2 regions: 1) Dominikia bernensis (>2000 reads); 2) Glomus indicum, Rhizophagus intraradices and Entrophospora infrequens (>1000 reads for each species); 3) Ambispora gerdemannii, Rhizophagus irregularis and Paraglomus laccatum (>500 reads for each species). While the culculation of the number of OTUs for the ITS1+ITS2 regions indicated that the most common species were as follows: 1) Rhizophagus intraradices (51 OTUs); 2) Paraglomus laccatum (37 OTUs); 3) Rhizophagus irregularis (33 OTUs); 4) Dominikia bernensis (32 OTUs); 5) Claroideoglomus claroideum (21 OTUs); 6) Funneliformis mosseae (17 OTUs), and 7) Acaulospora paulinae (16 OTUs). The most common genera were: 1) Rhizophagus and Dominikia (>3000 reads for >100 OTUs for each genus); 2) Glomus (>2000 reads for >50 OTUs); 3) Paraglomus and Claroideoglomus (>500 reads for >30 OTUs). The list of endemic species is presented in Table S8. The results showed that among the known species, 7 specific (endemic) species were found for the biotopes of the subalpine meadow, 4 – in the forest, and 17 – for disturbed ecosystems, pastures in the biotopes of the river valley.

For a more accurate comparison of biodiversity in the studied biotopes, the Shannon and Simpson indices were used (Table S9). The AMF species diversity indices calculated from the data obtained using two various markers differ from each other. All Hill numbers were higher for the data obtained for ITS1 than for ITS2. This was typical for both observed and extrapolated indices. The lowest observed diversity indices were noted for the forest ecosystem, the highest ones were for the river valley. The Hill indices (estimated by sample size extrapolation) also revealed the low diversity for the forest ecosystem. The extrapolated value of species richness was maximal for the meadow, but the confidence interval value was also very high for the meadow, so it was complicated to conclude which of the two biotopes had higher index. The extrapolated Shannon and Simpson indices were significantly (P<0.01) higher for the river valley and significantly (P<0.01) lower for the forest according to ITS2 data, but close to each other in all three biotopes according to ITS1 data. A comparison of the observed and extrapolated indices showed that for the river valley the identification of more than 96% of the species diversity was performed in this study with both markers.

The obtained results indicated that the analyzed Caucasus biotopes had different biodiversity and included a different number of specific species of AM fungi. Therefore, it is of an interest to analyze the correlations between biodiversity for individual STPs and agrochemical parameters of sampled rhizospheric soils.

3.4. Correlation Analysis of the Interlinks between AMF Species Richness with STP Height and Agrochemical Parameters of Sampled Rhizosphere Soils

The assumption of this study was the existence of a correlation of species richness with soil agrochemical indicators and STP height (altitude). The analysis was performed with application of rho, Spearman correlation (Figures S6 and S7). The results indicated the existence of a direct positive correlations of AMF biodiversity only with the pH value of soil acidity. The highest greatest correlation was found for OTUs in the ITS1 region (Figure S6). The absence of correlations other soil parameters, as well as with the altitude above sea level was discovered. Linear correlation coefficients (r) between AMF biodiversity parameters, plant biodiversity and agrochemical indicators of STPs are presented in Table S10. A proved (P<0.05) positive correlations were found between the indicators of AMF diversity with “Percentage of annual plants”, with “pHKCl” and in some cases with “Pi, mg/kg” (inorganic phosphorus available for plant nutrition).

4. Discussion

4.1. Applicability of Molecular Genetic Identification of AMF in the Study of Their Biodiversity

Our study focuses on to define AMF biodiversity in the North Caucasus, a temperate biodiversity hotspot for plants and animals, to detect new AMF hotspot. The correct identification of individual AMF species is possible by the morphological method. But AMF do not grow on culture media, so the accuracy of estimating the AMF species richness by this method is low. In this regard, we have been carried out identification using molecular genetic methods based on using Illumina MiSeq. The results are revealed that significant differences among the three analyzed biotopes for Glomeromycetes class, which included all known AMF species. An assessment of the suitability and relevance of the approach used for molecular genetic identification of AMF, its advantages and disadvantages will be discussed further on. We proposed new methodology to solve the problem of high polymorphism of AM fungi for ITS regions by constructing trees by genus and accompanied with manually-verified alignment (Figure S1). Because of the impossibility to construct a phylogenetic tree for the whole Glomeromycetes class for ITS regions, the most studies of AM fungi are provided with conservative genes (for example, the LSU region; [77]). The LSU region is also not fully suitable for constructing a complete reference tree for Glomeromycetes. In the investigation of Delavaux et al., this tree included only 174 sequences of AM fungi assigned to only 112 species [77], despite the fact that 345 species of AMF are known [30]. Thus, this reference tree is not enough to identify most of the species (more than ⅔ of known AMF species). Therefore, the assessment of AMF biodiversity and their identification only by LSU at the species level is not yet possible, and most studies are currently associated with the definition of AMF taxa to the genus level [14,22,17,77]. The SSU region is believed to be slowly evolving, and therefore not sufficiently variable to adequately identify AMF species [78–80]. Suggested in this study algorithm (Figure S1) allows one to eliminate the disadvantage of using polymorphic ITS regions, to carry out correct clustering for AMF identification at the species level. Currently, ITS1 and ITS2 are, as a rule, analyzed unequally [31]. Thus, in the NCBI GenBank database, individual sequences of the ITS2 region are the most frequent, and sequences of the ITS1 region are represented as part of the entire SSU-ITS1-5.8-ITS2-LSU region. According to the data of our study ITS1 analysis demonstrated a significant number of reads, and it can be assumed as a second barcoding region along with the ITS2 region. In all likelihood, the efficiency of primer annealing for the ITS1 and ITS2 regions differs, which is confirmed by the different composition of the identified AMF species; Tables S4 and S6). The assumption is that the analysis of ITS1 and ITS2 regions should complement each other, as if they were different genetic markers.

In the biotopes of the North Caucasus, for ITS1 and ITS2, the total number of AMF species was 50 (the number of clades of the species level with accurate identification) from 20 genera. Currently, 345 AMF species from 44 genera are known [30]. Thus, at least 15% of the world's AMF species diversity has been identified in the studied region of the North Caucasus (from 45% of genera). AM has been actively investigated for several decades, including in the North Caucasus and other regions of Russia. The plant communities of the Teberdinsky National Park, including arbuscular mycorrhizae, have been studied for many years by Onipchenko et al. [81–83]. Authors found out that about 74-77% of the higher plants of the National Park interact with AMF, forming a symbiosis, arbuscular mycorrhiza [83]. However, a detailed assessment of the diversity of AMF species in the North Caucasus has not yet been performed. Evaluation of biodiversity and comparison of different floral complexes (river valley, mountain forest and subalpine meadow) was a difficult task without the application of molecular genetic methods. The main approaches were at the level of classical floristics, geobotany and mycology known as methods of comparative morphological analysis based on species identification only by the morphology of rhizospheric extraradical AMF spores, but not by intraradical AM structures: vesicles, arbuscules and mycelium [84]. Meanwhile AMF peculiarities are concerned to their obligate status in relation to the host plant and fundamental impossibility to grow AMF spores in soils without plants. Therefore, it is very difficult to assess AMF biodiversity by the morphology of spores. On the other hand, the difficulties of molecular genetic identification of AMF by rDNA regions linked to the fact that the concepts of “organism” or “species” are not fully applicable to AMF [38], because AMF contains nuclei of different origin in their hyphae and spores (heterokariosis), i.e. essentially one fungus contains a set of polymorphic genotypes, and since the analysis of genetic diversity is evaluated by the diversity of rDNA, it more correct a definition that one AM fungi has a set of different ribotypes. It is appropriate to consider the community of nuclei within the arbuscular mycorrhiza in the studied STP as a peculiar pangenome, rapidly changing due to horizontal gene transfer [85,87]. It is also known that AMF are able to form anastomoses [88] and, accordingly, the exchange of genetic material. All this is the reason that different OTUs belong to the same AMF species according to the sequencing results, obtained in this and earlier studies [62]. Thus molecular genetic methods including NGS sequencing using Illumina MiSeq have advantages in assessing AMF biodiversity up to the species level in comparison with morphological methods of AMF identification (Table S11). Nevertheless only morphological methods provide reliable information for accurate species identification. And some disadvantages of NGS sequencing should be mentioned (Table S11). The main problem is the high genetic polymorphism of AM fungi.

The genetic mechanisms of heterogeneous nuclei interaction within a single pangenome have not been studied. We assume that the genetic diversity of OTUs revealed in this study with Illumina MiSeq application may be caused by three different phenomena. Firstly, it is a consequence of the heterogeneity of nuclei in hyphae and AMF spores. For example, one spore can contain up to 35,000 nuclei [91]. The reason is that АМ fungi (unlike other eukaryotes) lack the genetic bottleneck of a single-nucleus stage [92].

Secondly, the detected heterogeneity of the AMF genetic material may be a result of conflict of genomes, a “genomic shock” (known from studies on hybrids and allopolyploids of plants). At the rDNA level it is expressed as “nucleolar dominance” known as silencing of the rDNA of one of the parents, a high rate of accumulation of substitutions in the repressed genome and a gradual decrease in the proportion of rDNA of the repressed subgenome in the polyploid genome [93–95]. Such a phenomenon of nucleoli dominance is also found in AM fungi [96]. It is displaied by the appearance of a high number of singletons – unique and quasi-unique rDNA variants represented among reads by a minimum number of copies [93], similarly to what was found in this study. The provided identification of a larger number of OTUs for ITS1 due to the greater ITS1 variability, and, secondly, since AMF contain many nuclei that have different variants of both ITS1 and ITS2. Consequently, the same mycelium of the same AMF species may contain more variants of ITS1 than ITS2. The different variability of the regions (the rate of changes in ITS2 is less than in ITS1) and the different specificity of primers to them indicate that the complex analysis of the both regions significantly complements each other. Singletons oftenly contain chimeric sequences, but, its exclusion can lead to a significant loss of truly existing fungi species. This was estimated in the analysis by Baldrian et al., since 2.5% of 10 million singletons had a more than 97% similarity to fungi species [25].

The third mechanism that can lead to the increase in the diversity of OTU in the studied STPs (in addition to the real taxonomic diversity of the mycota of each site) is a parasexual process that naturally directs to the exchange of genes between nuclei during karyokinesis [97]. The common opinion was that mycorrhizal fungi are asexual. However, the preservation in their genomes of 85% of gene complexes associated with crossing over of sexual forms [98] fully admits the assumption that a hidden parasexual process is possible [84,99].

4.2. Comparative Analysis of AMF Diversity in Different Biotopes

As a result of the OTUs, as well as phylogenetic trees analysis, the AMF species were identified (Tables S4 and S6). The ratio of OTUs and AMF species in three biotopes (subalpine meadow, forest, river valley) is represented in Figure 5 and Figure 6 for the ITS1 and ITS2 regions, respectively. The largest number of specific AMF species was found in the river valley, and the least one was in forest biotopes for both ITS regions, including calculation of specific species without VT (Table S8). The highest number of species has been identified in all biotopes for the ITS1 region in comparison with ITS2.

Among the biotopes of meadows, forests and river valleys, both endemic and common AMF species for different biotopes have been identified (Table S8). The dominant AMF species in the North Caucasus ecosystems are (in descending order of the number of OTUs, species with the number of identified OTUs > 350): R. intraradices, R. irregularis, Entrophospora infrequens, Dominikia indica. Whereas the conclusions about the distribution of a taxon based on the number of OTUs are somewhat risky, it should be noted that they may be correct in the case of dominant and subdominant species. For example, such AMF genera as Rhizophagus and Dominikia are recognized as the most common ones (by the number of OTUs and reads), so it should be considered as cosmopolitan genera. According to literature data [100], R. irregularis and R. intraradices species have been identified in many ecosystems of the world. The Entrophospora genus should be mentioned separately. This is a monophyletic group which has recently been included in AMF [101], and later assigned to the new order of the Entrophosporales (according to some sources, along with Claroideoglomus; [39]). The study of the Paraglomus genus is of a particular interest because its sequences in the NCBI database for ITS regions are not enough for accurate species identification (there are only 3 species out of 9 known by morphology of spores in the database). All AMF of these genera have a significant number of reads in the data sets obtained for the STPs of the North Caucasus (Tables S4 and S6). Furthermore the list of the most represented genera in the ecosystems of the Andes [57] differs from genera in the biotopes of the North Caucasus (Rhizophagus, Dominikia, and Glomus; see Tables S5 and S7). The following genera can be found in the Andes: Acaulospora (>50% reads; the genus is characteristic of high zones), Claroideoglomus, Cetraspora, Rhizophagus [57]. In the Southern Central zone of Chile, the main AMF genera were Acaulospora and Glomus [42]. In Yucatan, Mexico, the main AMF genera associated with tree roots were the others: Glomus, Sclerocystis, Rhizophagus, Redeckera, and Diversispora [17]. In Mt. Taibai of Qinling Mountain, China, the main AMF genera were: Glomus, Septoglomus, Acaulospora, Ambispora, and Rhizophagus [102]. We can assume that mountain biotopes differ significantly in the composition of AMF genera.

Comparative analysis of biotopes revealed common features in the species composition for ecosystems of the same type (Table S4 and S6). In addition to the known species, we identified 52 species-level VTs assigned to 12 genera (for the ITS1 and ITS2 regions). Virtual taxa found in the ecosystems of the North Caucasus, a temperate biodiversity hotspot, are also important because it may represent new, previously unknown uncharacterised AMF species. So in the study in Yucatan, Mexico, Lara-Pérez et al. [17] identified 36 VTs, belonging to 9 genera (for 18S region) were identified. In the Serengeti National Park, Tanzania, Stevens et al. [21] identified 39 VTs (for the SSU region); in the High Arctic, the Zackenberg valley, Northeast Greenland, Rasmussen et al. [43] identified 29 VTs (for SSU region); in the Middle Caquetá River region and the municipality of Leticia, the Colombian Amazon region, Peña-Venegas et al. [20] identified 126 VTs (for SSU region), but the application of SSU as a barcode did not allow the identification of AMF to the species level. The composition of morphological and molecular genetic approaches can significantly expand the possibilities of AMF identification. For examples, in the Iron triangle region, Brazil, Vieira et al. [27] found as many as 62 AMF species from 18 genera (for region 25S).

According to experimental results, for any type of data calculation (reads, OTUs, the number of precisely defined AMF species, the total number of species taking into account the number of virtual AMF taxa at the species level, the values of biodiversity indices), forest ecosystems have the least AMF biodiversity, and river valley ecosystems, on the contrary, have the highest AMF biodiversity. However, according to the estimated species richness (Table S9), the subalpine meadow had equally high diversity with that for river valley meadows. Based on the analysis of over 2 million valid sequences with an average sequence length of 216 bp for each sample, it was shown that the subalpine meadow (2,635 m above sea level) on the western slope of Helan Mountain in Alxa Left Banner, Inner Mongolia Autonomous Region, had the highest diversity according to the Shannon index, the ACE and Chao1 indices, though it manifested the smallest Simpson index [49].

4.3. Reasons for Higher AMF Diversity in River Valley Biotopes

To elucidate the reasons for the taxonomic diversity of AMF, we analyzed a number of factors that could potentially be important. Such abiotic and biotic factors include: climatic conditions, agrochemical composition of soils, human activity, and composition of plant communities [46–49]. Provided multivariate analysis of variance revealed that the STPs of one biotope (subalpine meadow /forest /river valley) are clustered by the number of OTUs in the space obtained using multidimensional scaling (MDS-Spearman). Correlation analysis indicated a significant (P<0.01 for ITS1) positive relationship between soil acidity (pH_KCl) and AMF species richness (Figures S6 and S7), and in some cases of AMF biodiversity indicators with the level of available P for plants (“Pi, mg/kg”; Table S10). The effect of other agrochemical soil parameters, and also the altitude of the STPs was weaker or absent (Table S10).

The list of possible reasons for high AMF diversity in a river valley:

(1)

Human activity (livestock grazing, etc.) is much more developed in river valleys than in mountain forests and subalpine meadows [12]. The invasion of AM fungi by transferring spores on human shoes and hooves of livestock (cows, horses, etc.) is quite common [48,103]. However, some evidances indicate that ungulate grazing may be associated with the decrease of AM fungal abundance in soil [21]. Along with that AMF spores can migrate with water flows from the mountains to the river valley during erosion [103–105] and after sedimentation it actively enter into nonspecific symbiotic relationship. It is developed despite the fact that the organic reserves in these ecosystems are much higher because of livestock grazing, moreover the ecosystems themselves may have signs of soil degradation and are considered to be disturbed. Nevertheless, AM fungi spores, as a rule, are significantly larger (>40 microns) than the spores of many other fungi, so their distribution distance is relatively short [106]. According to Guo et al. [107], terrain slope can also affect AMF diversity. The biotopes of the river valley analyzed in our study were characterized by much more gentle slope in contrast to the biotopes of subalpine meadows and forests (Table S2). Therefore, it can be assumed that a flat slope will have a positive correlation with AMF biodiversity.

(2)

In conditions of intensive percolation water regime and good drainage, there is no stagnation of water and oxygen deficiency in soil, negatively affecting the development of AM fungi [108,109]. The percolation water regime can reduce the content of available phosphate (Pi) in the soil, which makes root mycorrhization an important adaptation for P uptake. It is the water regime is typical for river valleys in the North Caucasus. The correlation between the parameters of AMF biodiversity and “Pi, mg/kg” (inorganic P available for plant nutrition) in many cases was significantly (P<0.05) positive (Table S10). Our data is consistent with results of Guo et al. in terms of that the diversity of AMF had a positive correlation with available P content [107]. Nevertheless, large-scale studies have shown that AMF diversity and abundance decrease with phosphorus available in the soil [Ma et al., 2023]. This issue requires further consideration.

(3)

STPs altitude may be a factor influencing AM development. The analyzed biotopes of the forest and subalpine meadow located above the biotopes of the river valley were characterized by a reduced number of identified AMF at the species level (Figure 5 and Figure 6), although linear reliable correlations for the North Caucasus STPs were not found (Table S10). Our results are consistent with literature data [10,113–116]. Thus, AMF biodiversity was practically independent of altitude [50,54–57], or the correlations were negative [10,113–116]. However, this rule is not applicable for the Zackenberg valley in the High Arctic [43]: with an increase in altitude from 33 to 479 m (a small altitude above the sea level), an increase in AM fungi occurrence was observed. Perhaps the reason for the lack of correlation between altitude and AMF biodiversity in the North Caucasus is that there are no biotopes of alpine meadows in the analyzed STPs, which are characterized by lower temperatures. However temperature is considered an important factor for the development of AM [109,117,118]. The optimal average air temperature of the warmest month is +19 °C, but mycorrhizal colonization can be intensively increased with the frosty period not less than 2 months [119].

(4)

The biodiversity of AMF and the development of AM are affected by such factors as pollution, salinity, drought, extreme temperatures, CO₂, liming, acidity, etc. [120], as well as soil composition, altitude, species composition of plant communities, climatic factors [46–49]. But it is pH, along with temperature, that are the main factors determining AMF biodiversity [59]. pH can have an important direct effect on the growth and productivity of the AM fungus [Wang et al., 1993; Coughlan et al., 2000]. The positive correlation between AMF diversity and pH is mentioned in [107,121]. This is consistent with our data on a significant positive correlation of the pH_KCl with a number of biodiversity indicators (“Glomeromycota total reads”, “OTUs number for total species”, “total species number”, both for ITS1 and ITS2; Table S10; Figures S6 and S7).

(5)

The phenotypic diversity of OTU AMF is supposed to be under the effect of the phenotypic diversity of plants, which decreases with altitude in the mountains [122]. The colonization of root by AM fungi is not species–specific. The roots of one plant can be colonized by several species of AM fungi, and one AM fungus can colonize different plant species [123]. However, despite the absence of direct correlations between the diversity of AMF and the total number of herbaceous plants (Table S10), it might be expected a decrease in the spectrum of potential partners in the mutualistic symbiotic system, and thus affect AMF species diversity [124]. At the same time, annual plant species have a higher diversity of AMF than perennial plant species, and half of the currently identified AMF species may be more specific to one plant species [53]. Moreover, S. Horn et al. [50] demonstrated that the influence of biotic factors (interaction of AMF with plants) is more significant in comparison with the effect of abiotic factors on the composition of AMF genera. It is known that in the process of succession with an increase in the proportion of woody plants the density of AM fungi spores decreased [51,52]. Thus, higher abundance of annual plants (see “Percentage of annual plants” in Table S2) in river valleys in comparison with the biotopes of the subalpine meadow may be a key factor positively affecting AMF taxonomic diversity. Meanwhile, our studies confirm the relationship between the proportion of annual plant forms and the diversity of AMF. For instance, it was shown that the linear correlation coefficients were reliable (P<0.05; Table S10). The correlation coefficient between “Percentage of annual plants” and “Glomeromycota total reads” r = 0.76 and 0.81 (for ITS1 and ITS2, respectively), and the correlation coefficient between “Percentage of annual plants” and “OTUs number (for total species)” r = 0.67 и 0.77 (for ITS1 and ITS2, respectively). A weak correlation between the proportion of annual plants and AMF species diversity was also shown (r = 0.40 and 0.56 for ITS1 and ITS2, respectively; see Table S10). Similar results were obtained in [42]: the “natural grassland” ecosystem had the highest AMF species diversity among 20 ecosystems of interest. The opposite is also possible. It is shown that in the Teberdinsky National Park, experimental suppression of AM symbiosis always is followed by a decrease in the species richness and number of plants [125]. The observations of OTU diversity in the studied territories in a different season may provide new information, since the species composition and numerical ratios of different OTUs may alter due to the season [126,127]. Changes in the mycorrhization of plants by AM fungi thrughout a year in the Teberdinsky National Park have already been studied earlier [81], but their biodiversity has not been assessed.

4.4. Practical Application of the Results of AMF Biodiversity Research

Rhizophagus irregularis, Funneliformis mosseae and Gigaspora margarita are the species that are most often used in biological preparations to enhance plant growth [44,128–132]. Meanwhile, the results of this study showed that among 10 million of sequences the ones with F. mosseae and Gi. margarita are extremely rare, which requires special attention of researchers, because these fungi are isolated from soils much more often by morphological methods. The isolation of AMF from the analyzed STPs, we collected a large number of R. irregularis spores, a few F. mosseae spores, but failed to isolate Gi. margarita spores from the roots of plants in the North Caucasus (unpublished data). Universal primers are, perhaps, less suitable for sequencing these two species. Other authors reported that these species may be absent [22,23], or rarely occur in the sequencing results, despite the fact that they are identified by morphological methods [9,17,21], or not all species of Gigaspora and Funneliformis genera are found [16,20,27]. Thus, we can assume that the application of universal primers (for ITS1 and ITS2 regions) make it possible to assess AMF biodiversity. In the case of studying a single new species, it may be useful to analyze, for example, the LSU region [77] as another barcoding region.

5. Conclusions

The investigation of the North Caucasus region, in the Teberdinsky National Park and adjacent territories of Karachay-Cherkessia, revealed a significant number of AMF species were identified for the first time: 50 species from 20 genera in three biotopes, like subalpine meadow, forest and river valley. It represented at least 15% of the world species diversity of AMF (from 45% of known genera). Thus the biotopes of the North Caucasus should be considered as a native temperate biodiversity hotspot not only of plants and animals, but also of AM fungi. Modification of the molecular genetic identification of AMF showed the effectiveness of the proposed algorithm (Figure S1), including: 1) using Illumina MiSeq simultaneously for two barcodes (ITS1 and ITS2 regions); 2) inclusion of analysis for genera phylogenetic trees as a separate step as a part of species identification algorithm; 3) selection of sequences from the MaarjAM database as OTUs centroids, which made it possible to use about 15% of reliable singletons and thereby minimize losses of true AMF species. The list of species identified in all biotopes includes: Acaulospora paulinae, Claroideoglomus claroideum, C. lamellosum, Dominikia indica, Glomus hoi, Paraglomus laccatum, Rhizophagus irregularis, R. melanus, R. intraradices. The dominant species, R. intraradices, R. irregularis and D. indica The species specific (endemic) for certain biotopes were identified: 7 species – for biotopes of the subalpine meadow, 4 – in the forest and 17 – for disturbed ecosystems, pastures in the river valley biotopes (Table S8). 52 virtual taxa of the species level were also found, assigned to 12 genera. These taxa may include previously unknown AMF species. Correlation analysis revealed that the main reasons for the high biodiversity in river valleys can be as follows: 1) a higher proportion of annual plant species in river valleys in comparison with biotopes of subalpine meadows and forests, since it is annual plant species that have a higher diversity of AMF than perennial ones; 2) disturbance of river valley biotopes (active grazing on pastures) in contrast to undisturbed biotopes of forests and subalpine meadows, that is a sign of the process of regenerative succession, in which the role of AMF is high; 3) positive correlation of AMF diversity with pH in the soil. Our data prove the importance of future investigations in the North Caucasus region, in the Teberdinsky National Park and adjacent territories of Karachay-Cherkessia. Accumulated data will contribute to the expansion of our knowledge about the role of AMF involvement in the formation of natural biotopes.