code from "Unveiling the rich and functionally diverse Southern Brazilian Highland Grasslands soil funga for promoting conservation" v1.0

<div> <h3>contact: <a href="mailto:[email protected]">[email protected]</a></h3> </div> <h3>This is a static version of the same dataset available at <a href="https://github.com/LBMCF/SBHG_mycobiota">https://github.com/LBMCF/SBHG_mycobiota</a>.</h3> <p>All necessary code and data for replicating the results presented in the article entitled "Unveiling the rich and functionally diverse Southern Brazilian Highland Grasslands soil funga for promoting conservation" can be found in this repository. R and python scripts are named based on which analyses were perfomed, and throghouly documented for guiding the execution of each step in the analysis pipeline. </p> <p>Raw shotgun metagenomic data can be found under NCBI BioProject PRJNA1227467 along with additional environmental, geographic, and sampling metadata available within BioSamples associated with the aforementioned BioProject.</p> <p>From the methods section:</p> <div> <h3>2.1. Area characterization and soil sampling</h3> </div> <p>The SBHG sampled areas are located within the São Joaquim National Park in Southern Brazil, within the Campos Sulinos region (Supplementary figure 1a). The region's climate is classified as humid subtropical, without a marked dry season (Cfb, according to the Köppen-Geiger classification). The samples were collected between May 25 and 27, 2022, with an average temperature of 10.3ºC, with maximum of 14ºC, and minimum of 8ºC. The average annual precipitation ranges from 1600 to 1900 mm. The predominant soils in the study region are Lithic Humudepts and Lithic Udifolists (Soil survey staff, 2022). Soil was sampled in three sites with an average distance of 4.3 km from each other (Supplementary figure 1a). Within each of the three sites (S1, S2, and S3), twelve composite soil samples were taken, arranged in three transects separated equally by 15 m (Supplementary figure 1b). Each composite soil sample consisted of four simple soil samples collected in the edge of each sampling unit (SU), 1.2 m apart from the center. Adjacent soil samples were collected in each SU for physical-chemical analyses. A sterile garden shovel was used to remove ground vegetation from the sampling points and a sterile soil corer (⌀ 254 mm) was inserted 10 cm deep in the ground, focusing on the soil histic or umbric epipedon. Between each SU, disposable gloves were used, and all tools were sterilized with a 1% sodium hypochlorite solution and sterile water to avoid cross-contamination. Immediately after collection, the composite soil samples were stored in sterile plastic bags and acclimated at coolers with ice packs until storage at -20°C.</p> <div> <h3>2.2. Soil physicochemical analyses, abiotic and biotic variables</h3> </div> <p>Clay content was determined by the Bouyoucos method, pH value in water (1:1 soil–water ratio); SMP index, soil organic matter (SOM) content was determined using an auto-analyzer (LECO, TruSpec CHNS, St. Joseph, MI, USA), available P and K (extracted by Mehlich-1); and exchangeable Al, Ca, and Mg contents (extracted with 1.0 mol L−1 KCl), following Teixeira et al. (2017) and Tedesco et al. (1995). With the obtained data, the H + Al contents, potential cation exchange capacity (CECpH7.0), saturation of CECpH7.0 by Ca + Mg + K (V%), and saturation by Al (m%) were calculate according to CQFS-RS/SC (2016). Plant influence on fungal communities structure was operationalized as the measured green biomass volume (i.e. standing and/or prostrated live plants), calculated as the product of green biomass area and height, and plant species dominance (inverse Simpson index, Sommerfield et al., 2008) for each SU. As an indication for temperature stress, the heat load index (HLI) for each SU was calculated using terrain slope and aspect and latitude, following McCune & Keon (2002), considering an aspect of 315° would receive the highest diurnal heat load in the southern hemisphere. 2.3. Extraction, Quantification of metagenomic DNA, and Sequencing Each soil sample was vigorously homogenized by rubbing its content for ca. 5 minutes, and after removing all fine root material, 125-250 mg of soil were subsampled for metagenomic DNA extraction using the ZymoBIOMICS DNA Miniprep Kit (Zymo Research, California, United States), following the manufacturer's recommendations. The total DNA extractions were quantified using 2 µL of solution in a NanoVue absorption spectrophotometer (GE Healthcare, Chicago, United States). High-quality total DNA samples were sequenced in the Illumina HiSeq 2500 high-throughput sequencing platform (Illumina, California, United States), using a shotgun metagenomic approach with 150 bp long paired-end reads.</p> <div> <h3>2.4. Bioinformatics Analyses</h3> </div> <p>A PCR-free alternative to circumvent known biases in rRNA amplicon sequencing (Tedersoo et al., 2015; Tedersoo et al., 2022a; Rammitsu et al., 2021; Reynolds et al., 2022) was adopted, based on the miTags approach (Logares et al., 2014; Martínez-Pérez et al., 2022). Therefore, metagenomic sequence data was used to retrieve unbiased rRNA gene markers datasets for both 18S and ITS regions. Sequence quality (Phred score) was evaluated using FastQC v0.12.0 (Andrews, 2010) and MultiQC v1.13 (Ewels et al., 2016). Since the sequences showed high quality (average Phred score = 35.64, minimum = 33.94, and maximum = 36.45), no further quality filters were applied. All generated reads were filtered using SortMeRNA v4.3.6 (Kopylova et al., 2012) in pair-ended mode with UNITE 9.0 (Abarenkov et al., 2024) and SILVA SSU NR99 138.1 (Quast et al., 2012) as references, keeping reads with an e-value < 10-1. BLAST-like outputs from SortMeRNA were used to detect and cluster extracted reads into database derived OTUs (ddOTUs). Reads passing a query coverage and identity percentage quality filter of 90% with an alignment length bigger than 70 bp were maintained. Reads with more than one hit surpassing an e-value < 10-1 were assigned to the hit with smaller e-value and bigger percentage identity. Although we adopted a broad identity percentage to define ddOTUs, this accounts for variability between query and database subject as a result of ITS variability (Estensmo et al., 2021), with a 90% threshold mostly corresponding well to intraspecific ITS identity (White et al., 2013). Trophic modes and guilds were defined at the genus level using the FunGuild database (Nguyen et al., 2016). Only genera with probable and highly probable identities were maintained. Genera identified as belonging to more than one trophic mode or guild were considered for both. SqueezeMeta v1.6.3 (Tamames & Puente-Sánchez, 2019) pipeline was used for contig assembly, gene prediction and functional assignment. Trimmomatic (Bolger et al., 2014) was used to remove adapters and trim sequences, and contigs were assembled using Megahit (Li et al., 2016). Prodigal (Hyatt et al., 2010) was used to predict genes, and Diamond (Buchfink et al., 2015) was used for homology searches in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa & Goto, 2000). Gene read count was normalized and represented as transcripts per million reads (TPM). The SQMtools R package v1.6.3 (Puente-Sánchez et al., 2020) was used to retain only annotated reads belonging to Fungi.</p> <div> <h3>2.5. Statistical Analyses</h3> </div> <p>Relative abundance per sample was calculated as the division of ddOTUs read counts by sequence depth. Alpha diversity was then calculated within a interpolation-extrapolation framework at the block level using Hill numbers (Hsieh et al., 2016) with order q=0 (species richness), q=1 (effective Shannon index, or number of common species), and q=2 (effective Simpson diversity, or number of dominant species). As species richness and other alpha diversity measurements are only reliably comparable when assemblages are sampled at a similar effort (Chao et al., 2014), sample sizes were normalized at double the minimum sample size between all blocks. To estimate the effects of sampling area, edaphic variables, heat load, and grass species on fungal community structure, a PERMANOVA model was constructed with Bray-Curtis dissimilarity based on Hellinger-transformed ddOTU table, with a total of 9999 permutations. To reduce dimensionality and remove correlation among edaphic variables, the first principal component (PC1) from a Principal Component Analysis (PCA) was included in the PERMANOVA model. This axis explained 51.7% of the total variance in the measured edaphic data (Supplementary figure 2). PERMANOVA results are sensitive to multivariate dispersion [herein understood as variation in sample species composition within sites, (Anderson & Walsh, 2013)]; therefore, the model construction was followed by multivariate dispersion analysis. To visualize community dissimilarities between sites, ordinations were constructed using Non-Metric Multidimensional Scaling (NMDS) with Bray-Curtis dissimilarity, two dimensions (k = 2), and a maximum of 500 iterations. For all statistical analyses performed, a p-value < 0.05 was considered indicative of statistical significance. To check whether different sequencing depths between samples influenced PERMANOVA results and NMDS visualizations, the same aforementioned model construction and procedures were followed using a rarefied ddOTU table version. Data was rarefied using 1000 iterations at the sample level, at the lowest sequence count (484 sequences for ITS data, and 230 for 18s data). For each iteration, the resulting rarefied values were stored, and the average 1000 iteration values were used. All statistical analyses were