The communal lands of Nuevo Becal, Campeche, Mexico are located at 18°48’00” and 89°15’45” N 9°08’00” and 89°20’05” W and cover an extension of 52,800 ha (Fig. 1). The area is a continuum of the Calakmul Biosphere Reserve, the largest tropical reserve in the country, and an important biological corridor connecting in the south to the Maya Forest in northern Nicaragua and Belize. The area is a relatively flat relief, with an average elevation of 250 to 375 masl in the Zoh-Laguna plateau (Martínez and Galindo-Leal 2002). This region consists mainly of calcium carbonate platforms, and prevailing soils derived from the limestone parent material are shallow, young, and have poorly defined horizons (López-García et al. 1990). The predominant vegetation is represented by semi-evergreen forests that, according to tree height and phenology, are characterised by tropical tall-to medium-stature forests (≥ 25 m) with significant leaf loss between 25 and 50%, medium-stature forests with trees between 18 and 25 m high and seasonally inundated low-stature forests with trees between 8 and 15 m in height (Martínez and Galindo-Leal 2002; Henricus et al. 2007). The climate of the study area is characterised by average annual precipitation that ranges from 1100 to 1200 mm, with 600 mm of precipitation in the driest month, and an average annual temperature of 26°C, with wind from September to December (Martínez and Galindo-Leal, 2002). Common tree species include: Swietenia macrophylla, Metopium brownei, Pseudobombax ellipticum, Bursera simaruba, Dendropanax arboreus, Manilkara sapota, Lonchocarpus castilloi, Bucida buceras, Brosimun allicastrum,and Lysiloma bahamensis. Land use in the communal lands of Nuevo Becal includes forest management (selective logging), natural conserved forest, chewing-gum production, secondary vegetation ≥20 years, and agriculture. At present, 25,000 ha are under permanent forest management, and in 2016, this area was certified following the criteria and indicators for sustainable forest management according to the standards of the Forest Stewardship Council (FSC).
Description of the study sites
Forest management is characterised by selective logging based on diameter, with a minimum cutting diameter between 35 and 55 cm diameters at breast height (DBH) and cutting cycles of 25 years, low harvest intensity with 20–30% wood extraction in logging blocks of 500 h per year. The soil samples were collected in areas with tree stands six years after selective logging. The vegetation type is classified into medium sub-evergreen forest, and primary harvested species include: Swietenia macrophylla, Manilkara sapota,and Metopium brownie.Natural forests are characterised by the minimum management of gum harvests or conserved forests. Secondary vegetation is abandoned land lacking agriculture-dominant species for 20 yearsand a low-sub-perennial-flooded forest vegetation type.
Soil collection and physicochemical analysis
The sites were chosen according to the geographic delimitation of forest use, natural conserved forest, and secondary vegetation. Soil collection was carried out in November 2020 during the rainy season. We collected 25 soil samples in each environment, and randomly selected seven independent soil samples according to the best DNA quality and quantity to sequencing (n =21), and 14 samples per glomerospores isolation in each study area (n =42). To avoid edge effects, the samples were collected at least 300 m from roads. We collected soil samples with a trowel at a depth of 20 cm. Freshly fallen leaves were carefully removed after sample collection, and approximately one kilogram of soil was collected in polyethylene bags. Soil samples were homogenised, and with the aid of spatula, a 50 ml falcon tube was filled with soil for DNA extraction. Trowels and spatulas were thoroughly rinsed with H2O and cleaned with 70% alcohol among samples. The samples were kept refrigerated at −4°C and transported frozen to the laboratory, where they were stored in an −80°C ultrafreezer until processing.
The physicochemical characteristics of soil were determined using standardised methods. For organic matter and organic carbon, the methods reported by Walkley (1974) and Jackson (1976) were followed. Total nitrogen was evaluated following the methodology reported by Parsons et al. (1984), while calcium carbonate was determined following the protocol of Molnia (1974). The redox potential and pH were measured directly at the surface level of the sample, with the specific electrode for each parameter and an OAKTON mod. 300 series potentiometer.
DNA extraction and sequencing
DNA was extracted using the DNeasy PowerSoil Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions. The quantity and quality of the DNA were determined using 1.5% agarose gel and a Nanodrop 1000 (Thermo Scientific, USA). DNA was sent to the sequencing service of ZymoBIOMICS ™ to process and analyse the composition of bacteria (n = 21 samples) and fungi (n = 21 samples). The fungal 18S rRNA region was amplified using the fungus-specific primer pair AMV4.5NF (5’ AAGCTCGTAGTTGAATTTCG 3’) and AMDGR (5’ CCCAACTATCCCTATTAATCAT 3’) (Suzuki et al. 2021). For bacteria, the V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified. Briefly, the V3–V4 region of the 16S rRNA gene for bacteria and 18S rRNA for fungi was amplified using the Quick-16S ™ NGS Library kit (Zymo Research, Irvine, CA, USA) using the specific primers for each group. The library was quantified using real-time PCR in a mixture based on equal molarity to avoid the formation of chimeras. The final library mix was purified with the Select-a-Size DNA Clean and Concentrator ™ kit (Zymo Research, Irvine, CA) and then quantitated with TapeStation® and Qubit®. The final library was sequenced on Illumina® MiSeq ™ with Reagent Kit v3 (600 cycles). Representative sequences for each virtual taxa (VT) for fungi and amplicon sequence variant (ASVs) for bacteria are available on the NCBI website under Bioprojects PRJNA806153 and PRJNA808063, respectively.
To determine the genera or families of Glomeromycota taxa, phylogenetic analysis was performed. A representative sequence for each VT in the study and sequences from the MaarjAM database (Öpik et al. 2010), and NCBI (www.ncbi.nlm.nih.gov) were used for multiple alignment with MAFFT version 7.149bv (Katoh et al. 2002). The phylogenetic tree was assembled using the neighbouring method with the MEGA 7 programme (Tamura et al. 2011) and the Kimura-2 model (Kimura, 1980). The bootstrap method (Felsenstein 1985) was used with 1000 replicates to support the branches. Henningsomyces candidus (AF334916) sequences were obtained from the GenBank database (www.ncbi.nlm.nih.gov) and used as an outgroup. VT with ≤ 95% sequence similarity that did not form a cluster corresponding to a known genus were considered potential new species (Jiang et al. 2018).
Glomerospore isolation and morphological identification
The glomerospores were extracted from 50 g of dry soil by the sieving and wet decanting method with a sucrose gradient (Błaszkowski 2012). Glomerospores were grouped by morphology and mounted with polyvinyl alcohol-lacto-glycerol (PVLG) and PVLG + Melzer (v / v). Morphospecies were identified using the specialised literature (Schenck and Pérez 1990; Błaszkowski 2012). We classified them following Wijayawardene et al. (2020), and we accepted the generic term Rhizoglomus, proposed by Sieverding et al. (2014) instead of Rizophagus, proposed by Schüßler and Walker (2010). We also adopted the terms glomerospores and glomerocarps (Goto and Maia, 2006; Jobim et al. 2019). The slide-mounted reference samples were photographed with PVLG and Melzer with a Nikon 7200 camera coupled with a light microscope adapter (Leica, Wetzlar, Germany) with a micrometre to estimate spore size.
Bioinformatic and statistical analyses
To estimate the diversity of fungi and bacteria, as well as the sampling effort, species accumulation curves and coverage curves were created for each environment using the iNEXT Online programme. Alpha diversity (⍺-diversity: Chao1, Shannon, and Simpson index) and beta diversity (β-diversity) were calculated to assess the diversity of AMF and bacteria. To identify differences in alpha diversity and soil chemical properties, we used the ANOVA test with the R platform using version 4.2.0, implemented with the RStudio graphical interface (RStudio Team, 2015). One-way ANOSIM permutation tests were used to analyse differences in the communities of AMF and bacteria in conserved, managed, and secondary forests. Subsequently, we visualised the compositional data using non-metric multidimensional scaling ordination (NMDS). Multivariate analysis was based on Bray–Curtis distances of AMF and bacteria composition and abundance using 9999 permutations and the Bonferroni correction of probabilities. Similarity percentages analysis (SIMPER) based on AMF and bacterial abundance was used to identify the contribution of each species to differences between sites using PAST version 4.03 software (Hammer et al. 2001). The diversity and overlap of AMF virtual taxa, glomerospores, and bacteria in the conserved forest, managed forest, and secondary vegetation were visualised in a bipartite network with R package BIPARTITE using the number of reads or glomerospores per site as a proxy for interaction frequency (Dormann et al. 2008). To evaluate the potential associations between microorganisms and the physicochemical variables of the soil, a redundancy analysis (RDA) was carried out using R software version 3.5.2. Significance in all statistical analyses was established at α < 0.05.