Metabarcoding of Soil Fungal Communities from Perturbed Areas Reveal Drive Role of Environmental Factors in Microbial Diversity of Colombian Andosols


 Changes in soil fungal community caused by land use have not been sufficiently studied in South-American Andosols, considered globally as important food production areas. This study analyzed 26 soil samples of Andosols collected from locations devoted to conservation, agriculture and mining activities in the southeastern region of Antioquia, Colombia, to establish differences between fungal communities as indicators of the degree of soil perturbation. The study developed a novel heminested PCR with primers SSUmCf Mix, ITS4 and fITS7 to assess Arbuscular Mycorrhizal Fungi detection in a Illumina MiSeq metabarcoding on nuclear ribosomal ITS2 region. A non-metric multidimensional scaling allowed exploring driver factors of fungal community changes, while fitted Dirichlet-multinomial models and PERMANOVA tests allowed identifying the correlations between alpha diversity indexes and community dissimilarities, as well as the significance of land use effects on fungal community composition. Furthermore, response ratios were determined to assess effect size by land use over relevant taxa. Results suggest a good coverage of fungal diversity with a detection of 10,529 high-quality ITS2 sequences belonged to phylum Glomeromycota. The analysis shows strong correlations of Shannon and Fisher indexes with dissimilarities on fungal communities among land uses (r=0.94), related to variations in temperature, air humidity and organic matter contents that lead to significant responses in abundances of relevant orders (such as Wallemiales and Trichosporonales). The study highlights the rich fungal biodiversity of the tropical Andosols, their specific sensitivities to environmental perturbation factors, and the useful range of a metabarcoding approach to characterize soil fungal communities.


Introduction
Changes in the structure and composition of soil fungal communities are important indicators of soil quality as soil functions, related to nutrient cycling, ecosystem provisioning, and climate regulation, can decrease due to loss of fungal diversity (Brinkmann et  Tropical Andosols represent more than 50% of the global Andosol land area and are an important unit of global agricultural productivity due to their coarse texture, high porosity, low cohesion and tendency to accumulate more organic matter (Anda and Dahlgren 2020; Kassaye et al. 2020; Pushkareva et al. 2021). Speci cally, Andosols in central Colombia are considered complex microbial habitats and important water reservoirs, in view of their high spatial variation in phreatic zones, thickness of organic horizons and effective cation exchange capacity (Casamitjana et al. 2020). However, there is a lack of information to understand how human actions can affect the richness of fungal communities and their compositions in these Andosols. The analysis of soil fungal communities at local scale is the base of studies on their ecological behavior because local abiotic and biotic conditions, speci c to the bioclimatic region or habitat, affect soil fungal community composition (Alzarhani et al. 2019). Furthermore, assessing the effect of speci c land uses on individual soil microbial taxa is key to propose management practices of the soil microbial diversity to maintain or increase soil quality (Victorino et al. 2021).
Identi cation of fungal taxa associated with soils is necessary for understanding these complex communities (Pauvert et al. 2019). Techniques, such as metabarcoding based on sequences of the internal transcribed spacer (ITS) region, have been recognized as the universal barcode for a correct taxonomic identi cation of fungal community members with potential to detect rare species in samples of different nature, including poor sequenced taxa (Badotti et al. 2017;Forin et al. 2018). Characterization of fungal communities using metabarcoding of ITS region sequences have expanded the capability to identify land use effects that shape these communities across space and time, associated to factors such as leaf litter, soil nutrients, pH and organic matter content (Rosenfeld et al. 2018;Sommermann et al. 2018;Turley et al. 2020). Tropical soils are important units to study fungal ecological preferences and their driving factors, considering that diversity of fungal groups can peak in tropical ecosystems (Tedersoo et al. 2014). Previous studies tested laboratory protocols and obtained good ITS amplicon pools from soils typical of tropical ecosystems such as rainforest, dry forest, littorals, coasts and Andeans forest (Barnes et  Despite the extensive use of ITS metabarcoding to characterize fungal communities, particular settings on the primers used during PCR ampli cations are necessary to avoid universal primer bias and increase the identi cation of speci c fungal groups. For instance, ITS universal barcoding primers that are commonly used to characterize the whole fungal community of soil are considered suboptimal for Arbuscular Mycorrhizal Fungi (AMF) characterization (Berruti et al. 2017). However, the use of ITS speci c primers can improve the identi cation of AMF (Waud et al. 2014). Although SSU and ITS2 primers such as SSUmCf, fITS9 and fITS7 provide a good survey of AMF community structure and diversity patterns, molecular approaches developed for estimating soil AMF diversity are still weak at characterizing species and have discrepancies in the taxonomic composition (Krüger et al. 2009;Berruti et al. 2017). Since the combinations of fungal primers should be used to enhance identi cation of underrepresented groups by universal barcodes such as Glomeromycetes (George et al. 2019), evaluations of the impact of primer selection in AMF community analysis, in different environment samples, may be needed to draw more general conclusions regarding the primer selection and PCR protocol most suitable for metabarcoding studies on AMF communities (Suzuki et al. 2020). In addition, characterization of fungal communities in particular soil types by Illumina MiSeq metabarcoding and new ampli cation approaches may boost the molecular methods to identify AMF taxa as a result of new representative sequences obtained from AMF present in speci c host plants ( Barbosa et al. 2021).
This study characterizes fungal communities in tropical Andosols, developing a novel heminested PCR approach with primers SSUmCf Mix, ITS4 and fITS7 to assess AMF detection. We tested how features and relevant taxa of soil fungal communities vary among three land uses with different degrees of perturbation in the southeastern region of Antioquia, Colombia. Illumina MiSeq soil DNA metabarcoding on nuclear ribosomal ITS2 region was used to describe total fungal diversity and in particular those of AMF taxa of sites: (1) relatively undisturbed within conservation areas; (2) subject to pressure by agriculture activities; and (3) subject to considerable pressure by mining activities related to clay extraction.
The features of fungal communities compared in this study allows establishing the effect size of land use on soil biodiversity, based on the hypothesis that changes in features of fungal communities re ect the effect of different degrees of perturbation. In the frame of an increasing demand to understand the factors driving fertility and microbial diversity of soils, the results revealed how environmental factors and soil properties change among different land uses and alter the key fungal components.

Study site and characterization of land uses
The eld research took place during a rainy season in three municipalities (La Ceja, El Retiro and Rionegro), in the southeastern region of Antioquia, Colombia (Fig.1). The municipalities are located in the central Andes mountain range at altitudes between 2300 and 2700 m.a.s.l.. Pre-montane humid forest and mosaics of pasture-crops are typical ecosystems in the region, which exhibits a unimodal temperature cycle with an average temperature of 21°C and an annual bimodal precipitation cycle with an average monthly precipitation of 387 mm. Within each municipality, the study collected samples of Andosols topsoil in each of the three targeted land uses: unperturbed areas for conservation (UnPA), areas managed for agricultural activities (AgA), and areas of mining extraction (MeA).
A point-transect sampling was carried out along 100-m transects delineated from the centroid of the polygon that encompassed each land use. Three equidistant points (at 50 m) were selected along each transects. Five subsamples (soil cores at 20 cm depth) were taken within a 2 m radius of each point to make a composite soil sample, using a soil auger with a diameter of 5 cm. 27 composite soil samples were collected and stored in plastic bags at 4°C for further analyses. Environmental variables (soil temperature, air relative humidity, dew point temperature and barometric pressure) and soil variables (pH, total dissolved solids and electrical conductivity) were measured with a Kestrel 5500 Weather Meter and a Lamotte TRACER 1766, respectively. Organic carbon was quanti ed after by wet digestion method (Benbi 2018), and soil moisture following thermo-gravimetric technique (Susha Lekshmi et al. 2014).
DNA extraction and PCR ampli cation DNA was extracted, from 0.25 g of soil, previously sieved through 2mm pore size sieve, with DNeasy PowerSoil kit (Qiagen, Hilden, Germany), following the manufacturer's protocol. To improve the number of AMF sequences detected, the ITS2 nuclear ribosomal region was ampli ed using Invitrogen Platinum HotStart PCR Master Mix (Thermo Fisher Scienti c) by means of a heminested PCR approach. In the rst PCR, the entire ITS region (ITS1-5.8S-ITS2) was ampli ed with the fungal primers SSUmCf Mix (Krüger et al. 2009) and ITS4 (White et al. 1990). The cycling conditions were an initial step at 95°C for 5 min, 35 cycles at 95°C for 30 s, 54°C for 45 s, 72°C for 1 min, and a nal extension step of 72°C for 7 min.
The obtained PCR products (c.ca 750 bp) were checked on 1.2% agarose gel, diluted at 1:20-1:100 and used as template in the heminested PCR targeting the ITS2 region, with primers fITS7 (Ihrmark et al. 2012), and ITS4 endowed of the following Illumina overhang adapter sequences: forward overhang: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-and reverse overhang: 5′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-. The heminested PCR cycling conditions were as follows: an initial step at 95°C for 15 min, 27 cycles at 95°C for 30 s, 57°C for 30 s, 72°C for 30 s, and a nal extension step of 72°C for 7 min. PCR products were puri ed by means of Wizard SV Gel and PCR CleanUp System (Promega), quanti ed using Qubit 2.0 (Thermo Fisher Scienti c, Waltham, MA, USA) and sent to IGA Technology Services S.R.L (Udine, Italy) for Illumina MiSeq platform sequencing (2 x 300bp).

Bioinformatics and statistical analysis
Raw Sequence data were processed using the Quantitative Insights Into Microbial Ecology 2 software package (QIIME2 version 2020.6.0) (Bolyen et al. 2019). To avoid poor quality sequences, reads were truncated (> 280 bp for forward,> 265 bp for reverse reads). The classi er adopted to generate the taxonomic assignment was UNITE community 2020: UNITE QIIME release for Fungi version 20.02.2020. Sequences with ≥ 97 similarity were assigned to the same Operational Taxonomic Unit (OTU). The fungal sequencing data were uploaded to the NCBI SRA database under accession number PRJNA779046.
The generated dataset and metadata were used to create a phyloseq object with the R package "phyloseq" version 1. To identify taxa with a relevant role in the fungal communities, the number of samples in which each species occurs (expressed as average prevalence, aP) was estimated and the most frequent species in every land use were identi ed. A permutation-based test of Multivariate Homogeneity of group dispersions (PERMUTEST) followed by Permutational Multivariate Analysis of Variance (PERMANOVA) (each one with 999 permutations) was performed to assess the variation of environmental attributes, soil properties, fungal community composition (at different rank level) and alpha diversity indexes among land uses and municipalities. With PERMUTEST the assumption of equal data dispersion among samples was assessed.
Afterwards, abundances of taxa signi cantly different among land uses were highlighted by generating a heat-map diagram based on Bray-Curtis similarity coefficients (Pang et al. 2019). These coe cients were used to group soil samples with a hierarchical clustering analysis. Then, adopting a non-metric multidimensional scaling (NMDS) a gradient analysis was developed to highlight the driving role of environmental factors and soil properties on the ordination of those taxa abundances and alpha diversity indexes (Chao1, ACE, An unpaired two-sample Wilcoxon test was used to evaluate whether the Absolute Response Ratio (ARR) means were different between comparisons UnPA vs AgA and UnPA vs MeA.

Results
Sequencing and soil fungal community composition The PCR heminested approach produced a set of 3,524,761 raw sequences, with 353,312 high-quality ITS2 sequences, averaging 11,330 sequences per sample, retained after the denoising process. Twenty-six samples were used for downstream analyzes. Despite the higher number of valid sequences, land use MeA presented a lower number of OTUs. On the contrary, UnPA and AgA showed the highest number of OTUs (Table.1). Rarefaction curves suggest that a plateau was reached (Fig.2), indicating that the sampled environments have a good representation of the fungal community, with the identi cation of 489 OTUs, 9 phyla, 29 classes, 53 orders, 117 families, 158 genera and 174 species.
Among the 489 OTUs, 213 belonged to the Basidiomycota, 167 to the Ascomycota, 25 to the Mucoromycota, 17 to the Mortierellomycota, 16 to the Glomeromycota, and 41 to unassigned phyla. The heminested PCR approach allows identifying 10,529 (2.98%), high-quality ITS2 sequences of phylum Glomeromycota (4 orders and 6 families). Claroideoglomus (Family Claroideoglomeraceae) was the unique genus assigned. UnPA (5,302 sequences) has the highest number of Glomeromycota sequences, followed by MeA (2,709 sequences) and AgA (2,518 sequences) land uses. Some families such as Ambisporaceae and Acaulosporaceae were only present in the mining sites of La Ceja (LC-QC) ( Table.2).
UnPA contained 116 species, AgA 68 species and MeA 47 species. The most frequent species in every land use have similar prevalence (Table.3), with the frequent species in MeA land use being the less prevalent. Saitozyma sp. was the most frequent species in UnPA, whereas Trichoderma asperellum and Agrocybe pediades were the most frequent in AgA and MeA land uses. Malassezia sp. was prevalent and highly abundant in both perturbed land uses.

Soil fungal communities among land uses
The municipalities did not show signi cant differences at phylum and order rank. However, fungal taxa varied more strongly with land uses at the rank of order followed by the genera. Contrasting fungal genera abundances among land uses were only signi cant

Fungal diversity responses to variations in soil properties and environmental factors
The studied Andosols have an acidic to weakly acidic pH range (from 5.3 to 5.7), with some parameters such as organic carbon content and soil moisture decreasing with the degree of perturbation. On the contrary, other parameters such as electrical conductivity and total dissolved solids increase with the degree of perturbation. Similarly, environmental parameters like relative air humidity and barometric pressure decrease in perturbed areas and other parameters such as soil temperature increase in those areas (Table.4). In AgA samples, pH and soil temperature decreased. In general, soil properties and environmental variables show a signi cant variation among land uses (R 2 =0.3182, p=0.0002).
The NMDS dissimilarity distance matrix with Bray-Curtis, performed to assess variability in community composition among land uses and speci c activities, revealed that fungal communities of MeA were distinct from the other two land uses, observing a closer similarity between fungal communities from UnPA and AgA. The communities order composition positively correlates with relative air humidity and organic matter and negatively correlates with soil temperature and dew point temperature (Fig.6a). Three groups in the sample ordination, based on alpha diversity indexes, showed positive correlation with organic matter and negative correlation with pH ( Fig.6b). The subset of indexes with the maximum rank correlation with samples dissimilarities was shaped by Shannon and Fisher index (r=0.94). Each vector had a substantial effect on the sample gradient ordination (p=0.001). However, every subset of indexes strongly correlates (r=>0.70) with dissimilarities in the communities.
Effect size measure of land use UnPA and MeA showed a signi cantly high contrast in the ARR or effect size of differences observed in fungal community features and parameters assessed (Table.5). In general, land use intensity negatively affected abundances of prevalent fungal species, dominant order, alpha diversity index and environmental variables. The factors signi cantly affected by MeA and AgA land uses were the abundances of taxa belonging to order Wallemiales, Trichosporonales and organic matter. Despite this, the RRi soil moisture was signi cant in contrasting UnPA and MeA land uses, but not signi cant in contrasting UnPA and AgA land uses (Fig.7).

Soil fungal community in Colombian Andosols
The average of high-quality ITS2 sequences and the percentage of sequences retention obtained with the heminested PCR approach were comparable with previous sequencing efforts of fungal communities in similar tropical soils. In this study, Glomeromycota was a phylum rarely observed. However, among the found small set of high-quality ITS2 sequences of this phylum, 14 taxa were discriminated. Previous fungal community characterizations in comparable environments have reported less AMF diversity values. For example, Landinez-Torres et al. (2019) tested in Colombian Inceptisols (located at 2,900 m.a.s.l) an ITS1 region ampli cation by primers BITS and B58S3 identifying for Glomeromycota phylum (among 8 taxa) only 2 orders, 2 families (Glomeraceae and Acaulosporaceae) and 2 species (Glomus mosseae and Entrophospora infrequens). The low number of species of AMF discriminated through the heminested PCR approach tested can be ascribable to limits of universal primers targeting the ITS region. Nevertheless, nested PCR approaches could be preferred using mycorrhizal speci c primer for detection of AMF (Kumar 2018).
For instance, a similar PCR heminested approach with primers fITS9-ITS4 to amplify ITS2 region was useful to improved discrimination of AMF genus in a general fungal community characterization (Victorino et al. 2021).
In the Colombian Andosols analyzed in this study, sequences of phylum Basidiomycota were frequently present, especially in mining Furthermore, clustered samples from UnPA and AgA showed a high abundance of order dominant Hypocreales (phylum Ascomycota) and Agaricales (phylum Basidiomycota). In UnPA land use, factors such as higher air humidity, lower soil temperature, wild vegetation and the absence of pesticide applications can promote the abundance of the prevalent species belonging to Hypocreales order. For example, fungi belonging to Metarhizium genus are endophytic entomopathogenic whose preferences for habitat selection respond to temperate environments, vegetative cover, high humidity and associations with insect hosts (Guerrero-Guerra et al. 2013; McGuire and North eld, 2020). Likewise, Fusarium solani abundances are robustly associated with factors such as high soil moisture and presence of entomopathogenic nematodes (Wu et al. 2019). The higher prevalence of Saitozyma sp. (order Tremellales) in UnPA can relate to the litter observed in the sampled environments (Pre-montane wet forest and Forestry crops of Pinus sp.) as the Saitozyma genus is an important yeast able to incorporate carbon into soils. The high abundance of Saitozyma is typical of soils with abundant litter (Mašínová et al. 2017). In AgA land use, the high abundance of Ascomycota sequences was con rmed by the prevalence of Trichoderma asperellum, a fungal species with distribution patterns associated with typical Andean crops and coniferous forest Soil fungal alpha diversity was the best ordination model of samples gradient as a function of the degree of soil perturbation. The patterns obtained for every alpha diversity index among land uses indicated that fungal richness was higher in unperturbed areas and decreased conforming perturbation degree increased in each land use. This qualitative divergence between perturbation degree and fungal alpha diversity resulted in a signi cant association of dissimilarities in alpha diversity communities with the Shannon and Fisher index. Land use effect size over these alpha diversity indexes suggest that the primary fungal community features affected by land use relate to the ratio of the number of taxa and the abundance of those taxa (Fisher index ARR average= 0.88) followed by the taxa richness (Shannon index ARR average= 0.56). In the validation of those observations, the sequence abundances of the most relevant taxa (such as prevalent species and dominant order) were negatively affected by land use change, despite high average ITS valid sequences in agricultural (12,486 sequences) and mining areas (15,029 sequences). Furthermore, species richness downturned from 116 species identi ed in unperturbed areas, to 68 in agricultural areas and 47 species in mining areas.
The fungal community homogeneity at local scale observed in this study can represent a stable microbial consortium whose variability depends mainly on environmental and physic-chemical gradients in soil. The results support the hypothesis that relevant fungal taxa can be a useful tool for gaining insights into the structure of the soil microbiome, since relevant taxa in perturbed areas facilitate the identi cation of responses of soil microbiomes to environmental factors (Xiao et al. 2021). Likewise, the stronger response of alpha diversity indexes to the degree of perturbation found in this research is an important indicator of the fungal community perturbation patterns, rather than the abundance pattern of soil fungi taxa (Wang et al. 2015).

Conclusions
The species composition and structure of fungal communities in the assessed Colombian Andosols reveal the richness of soil fungal species compared to similar soil types. Our results provide a validation of the heminested approach using the primers set SSUmCf Mix-ITS4 and fITS7-ITS4 and Illumina MiSeq platform sequencing, to investigate fungal communities in Andosols. This approach has proven to be at least as accurate and comparable to other PCR based methods used for soil fungal community characterizations. In particular, the discrimination observed at rank of order and families indicate that the ampli cation of ITS2 region could be helpful to unravelling Glomeromycota diversity with major resolution.
The relevant taxa such as prevalent species and dominant order Hypocreales, Agaricales, Wallemiales and Trichosporonales were good indicators of the land use effect over soil fungal community. Factors such as air humidity, soil temperature and soil organic matter contents were more important as drivers than other edaphic factors that determined order composition dissimilarities among Data availability: The datasets generated and analyzed during the current study are available from the corresponding author upon publication. Order Archaeosporales  74  4707  237  568  1436  402  60  953  1470   F_Ambisporaceae  0  0  0  0  0  0  0  0  75   F_Archaeosporaceae  0  86  51  497  108  353  0  0  847   F_Unidenti ed  74  4621  186  71  1328  49  60  953  548 Order