Study area
The study is applied to an agricultural area on the southeastern edge of the city of Milan (
Figure 11), which has been the focus of recent strategic transformation plans aimed at establishing an agroecological experimentation laboratory,18 integrated into the Living Lab Milano Porta Verde system:19 an open-air laboratory dedicated to the socio-environmental regeneration of the area. The study area covers approximately 100 ha, part of the South Milan Agricultural Park - Vettabbia Valley system, an ancient strategic axis linking city and countryside.20 The area is predominantly composed of agricultural lands (arable crops, horticulture, permanent and fallow grasslands) and semi-natural areas (young reforestation). Since 2019, the implementation of 2 ha under productive agroforestry management - a multilayered agro-silvo-pastoral system inspired by regenerative practices - has begun (
Figure 1). The study area belongs to an alluvial context (Po Plain) dominated by Cambisols (WRB Classification System21), and Gleysols strips along water courses. A loamy-sandy texture predominates: loose and relatively permeable soils, easily workable, typical of alluvial contexts. Soil samples were collected 16–17 April 2023 for DNA analysis while samples for soil abiotic analysis were collected 27 June 2023. Four land use types were investigated as treatments: crop field (crop), uncultivated area (uncultivated), wooded area (reforested), agroforestry system (agroforestry) (a summary of management history (current management, previous land use, products) is included in
Table 1). Their selection was based on:
• the intention to compare the agroforestry system with different local land uses, yet still of agro-silvo-pastoral or semi-natural type;
• the opportunity to compare sites characterized by the same pedological substrate (same cartographic Unit22).
Table 1. Sampling sites main features.
| Land use type | Abbr. | Current management | Current management start date | Previous land use | Products |
|---|
| Annual crop field | Crop | Conservation agriculture | About 25 years ago | Conventional crop field | Corn (2–3 yrs) in rotation with winter wheat (1 yr) and soyabean (1 yr) |
| Wooded area | Reforested | Semi-natural
| 2018 (planting) | Uncultivated area | No products |
| Agroforestry system | Agroforestry | Regenerative agriculture | 2019 (planting) | Conventional crop field | Fruits, berries, eggs |
| Uncultivated area | Uncultivated | Mowings | Stable for decades | Crop field | Hay |
In each site, sampling was conducted at selected nodes of an 8 × 8 grid with a 20 m mesh (
Figure 11): we sampled randomly in each grid, for a total of 3 soil sampling points for each site, avoiding marginal points, which often coincided with the edges of the site (roads, ditches, tree rows) and were more likely to be subject to disturbances that could influence the analytical results. Points that coincided with specific sources of heterogeneity within the site were also excluded, such as those located very close (< 2 m) to tree rows crossing the uncultivated area. The sampling was carried out by removing litter and picking soil slices through a spade at a depth of 0–30 cm with fixed thickness.
Table 2 reports the soil parameters investigated. Abiotic parameters were measured at Labs & Technological Services AGQ, S.L. (https://agqlabs.com/en/).
Table 2. Abiotic parameters: complete names and abbreviations used in figures.
| Abiotic parameter | Abbreviation |
|---|
| Sand % | Sand % |
| Clay % | Clay % |
| Silt % | Silt % |
| Total Nitrogen mg/kg | Total N mg/kg |
| C/N Ratio | C/N Ratio |
| Organic Matter % | OM % |
| OC % | OC % |
| Available Phosphorus Olsen mg/kg | Avail. P Olsen mk/kg |
| pH | pH |
| Electrical Conductivity (μS cm −1, 20 °C) | EC μS/cm |
| Cation-Exchange Capacity meq 100 g −1 | CEC meq 100 g −1 |
| Base Saturation % | Base Saturation % |
| Exchangeable Calcium meq/100 g | Exch. Ca meq/100 g |
| Exchangeable Potassium meq/100 g | Exch. K meq/100 g |
| Exchangeable Magnesium meq/100 g | Exch. Mg meq/100 g |
| Zinc (DTPA) mg/kg | Zn (DTPA) mg/kg |
| Manganese (DPTA) mg/kg | Mn (DTPA) mg/kg |
| Iron (DTPA) mg/kg | Fe (DTPA) mg/kg |
| Boron (mg/kg) | B mg/kg |
| Copper (DTPA) mg/kg | Cu (DTPA) mg/kg |
| > 2 mm fraction | > 2 mm fraction |
Metagenomic sequencing library preparation, sequencing and analysis
The next generation sequencing experiments, which included quality control and the initial bioinformatics analysis, were executed by Genomix4life S.R.L. (Baronissi, Salerno, Italy). DNA quality and quantity were checked with the NanoDropOne spectrophotometer (Thermo Scientific, Waltham, MA, USA) and the Qubit Fluorometer 4.0 (Invitrogen Co., Carlsbad, CA, USA). 16S amplification was performed using primers targeting the hypervariable V3 and V4 region of the 16S rRNA gene 314F(5’-CCTACGGGNGGCWGCAG-3′) and 806R (5’-GACTACHVGGGTATCTAATCC-3′),23 while ITS amplification used the ITS3 – ITS4 primers (ITS3f 5′- GCATCGATGAAGAACGCAGC -3’ ITS4r: 5’-TCCTCCGCTTATTGATATGC -3′). PCR amplification and library preparation followed the Illumina 16S Metagenomic Sequencing Library Preparation protocol (Illumina, San Diego, CA, USA). Briefly, a two-step PCR protocol was used: the first PCR amplified the V3–V4 region of the bacterial 16S rRNA gene using locus-specific primers containing Illumina overhang adapters, while the second PCR attached dual indices and sequencing adapters using the Nextera XT Index Kit (Illumina). Amplifications were performed in duplicate for each sample to reduce PCR bias and the resulting products were pooled prior to library preparation. A negative control, including all reagents used in amplification and library preparation but excluding DNA template, was included to monitor potential contamination. The negative control included during library preparation did not generate a detectable library profile. During the sequencing run no reads were obtained from the negative control sample. Therefore, no sequences from the negative control were present or required filtering in the downstream analyses. Libraries were quantified using a Qubit fluorometer (Invitrogen Co., Carlsbad, CA, USA) and pooled to obtain an equimolar concentration of 4 nM for each index-tagged sample, with the addition of the PhiX Control Library. The pooled libraries underwent cluster generation and sequencing on the MiSeq platform (Illumina) using a 2X250 bp paired-end format. The amplicon sequencing data were analyzed using the 16S Metagenomics application available in the Illumina BaseSpace Sequence Hub. This pipeline performs read processing and taxonomic classification using a k-mer–based Naïve Bayesian classifier derived from the Ribosomal Database Project (RDP) algorithm.24 The RefSeq RDP 16S v3 database was used for taxonomic assignment of 16S rRNA sequences. For the ITS datasets, the same analytical framework was used; however, taxonomic assignment was performed using the UNITE reference database, which is specifically curated for fungal ITS sequences.25 Therefore, the pipeline itself was not fundamentally different, but the reference database used for taxonomic classification was adapted to account for the differences between bacterial 16S and fungal ITS markers. The UNITE Fungal ITS Database v7.2 is based on FASTA from https://doi.org/10.15156/BIO/587475 Includes singletons set as RefS (in dynamic files).
Statistical analysis
Raw taxa abundances were first processed by normalizing the values for each sample to relative abundances (ranging from 0 to 1). To visualise the composition of the most common taxa at genus-level, we generated heatmaps showing the relative abundance of the 25 most abundant fungal and bacterial genera across the 12 soil samples.
To quantitatively assess the composition of soil communities, we calculated two alpha-diversity metrics: observed richness (number of unique species) and Shannon diversity index.
These metrics were computed separately for fungal and bacterial taxa, based on species-level metagenomic assignments from each of the 12 soil samples.
Richness was defined as the count of distinct taxa detected per sample. The Shannon diversity index was calculated as:
(1)
where
N is the number of observed taxa in a sample and
p i is the relative abundance of the
ith taxon. To visualise the ecological indices, we created barplots of richness and Shannon’s index across all 12 samples and boxplots across the different land-uses (crop, reforested, agroforestry, uncultivated), for both fungal and bacterial species. The ecological indices were computed using custom Python (v3.12.4) code, while the plots were created using the seaborn (v0.13.2) and matplotlib (v3.8.4) libraries.
To complement diversity and abundance-based metrics, we also assessed the overlap of microbial taxa across land-use types using UpSet plots.26 These diagrams summarize taxa overlap among land uses: horizontal bars indicate the total taxa per land use, and vertical bars indicate the number of taxa in each intersection defined by the filled-dot combinations. Single-dot intersections represent land-use–unique taxa, whereas multi-dot intersections represent shared taxa. A taxon was considered present in a specific land-use type if it was found in one or more replicates. These plots visualize the intersections of species sets and were created using the UpSetPlot (v0.9.0) package in Python.
To assess how soil communities and environmental parameters vary across different land-use types, we performed comparative analyses focusing on diversity metrics, taxonomic composition, and abiotic factors. The 12 soil samples were grouped according to their land-use classification: crop, reforested, agroforestry, and uncultivated (3 replicates per group).
Richness and Shannon diversity values and relative abundances of fungal and bacterial species were compared across land-use types using boxplots. Statistical differences were assessed via one-way ANOVA followed by Tukey’s Honest Significant Difference (HSD) post hoc test.27 Significance was defined at p < 0.05. Given the large number of species detected, only those exhibiting statistically significant differences between land-use types were visualized. Assumptions for parametric testing were assessed by applying the Shapiro–Wilk test for normality and Levene’s test for homogeneity of variance across land-use types. We note that, given the small number of biological replicates (n = 3 per group), these diagnostics must be interpreted cautiously. To evaluate compositional turnover among samples and between land-use types, we computed beta diversity using the Whittaker approach. Starting from species-level relative abundance tables, we converted abundances to presence/absence by considering a taxon as present when its relative abundance exceeded a fixed threshold (p > p0). For each land-use type g, we computed the mean within-sample richness (alpha diversity) as
(2)
where
is the number of taxa present in replicate k of land use g and
is the number of replicates. We then computed the pooled richness (gamma diversity) within each land use as
(3)
where
denotes the set of taxa present in replicate k. Whittaker beta diversity within each land use was defined as
(4)
with larger values indicating greater turnover (i.e., lower overlap) among replicates within the same land-use type. To summarize turnover between land-use types g and h, we first defined the pooled sets
and
, then computed
,
, and
. Pairwise beta diversity between land uses was calculated as
(5)
where higher values indicate greater compositional turnover between pooled land-use communities.
Abiotic soil parameters were compared across land-use types using the same procedure. Electrical conductivity (EC) and boron (B) included left-censored observations (reported as <0.7 μS/cm and < 0.5 mg/kg, respectively). These values were set at one-half the respective detection limit for the analysis. Again, we represented only the parameters where significant differences were found. All visualizations and tests were implemented using custom Python scripts with seaborn (v0.13.2), scipy (1.12.0), and statannotations (v0.7.1).
First, we assessed samples similarity computing pairwise Pearson correlation matrices among the 12 soil samples using three different variable sets: bacterial species abundances, fungal species abundances, and abiotic parameters. These correlation matrices were visualized as heatmaps paired with hierarchical clustering.
Moreover, to identify and visualize the strongest associations among abiotic soil parameters, we computed pairwise Pearson correlation coefficients using the full set of measured variables for each one of the four land-uses types. To reduce visual clutter, we implemented a filtering strategy to retain only the most strongly correlated variable pairs.
Specifically, we first computed the Pearson correlation matrix across all abiotic parameters, using the absolute values of the coefficients. We then ranked all unique variable pairs by the magnitude of their correlation coefficients and selected the top 20 strongest associations (excluding self-correlations). The set of variables involved in these top correlations was extracted, and a reduced correlation matrix was constructed using only this subset. This filtered matrix was visualized as a heatmap with hierarchical clustering.
This analysis was carried out using custom Python (v3.12.4) scripts, while the heatmaps were created using seaborn (v0.13.2) package.
Before performing Multivariate Analysis, all variables were scaled to zero mean and unit variance using the StandardScaler function from the sklearn (v1.4.2) package. To reduce the dimensionality of the abiotic dataset and identify environmental gradients, we applied Principal Component Analysis (PCA) to the standardized abiotic parameters using the PCA function from sklearn (v1.4.2). The first two principal components were used to visualize the samples in reduced space, with samples colored by land-use type. Convex hulls were drawn around groups to highlight clustering by land use using the ConvexHulls function from scipy (v1.12.0).
To investigate the relationships between soil community composition and abiotic variables, we performed Redundancy Analysis (RDA) using species-level abundance data for fungi and bacteria as response matrices and a selected subset of abiotic parameters as explanatory variables. Because the number of abiotic variables exceeded the number of samples, we retained as explanatory variables only a selection of abiotic parameters. We first selected candidate predictors based on their contribution to the main abiotic gradient captured by PCA, retaining variables with an absolute loading >0.25 on the first principal component (PC1). To limit redundancy among predictors, we then assessed pairwise correlations among the selected abiotic variables and removed highly collinear variables using a correlation threshold of |r| = 0.9. The final set of abiotic predictors used in the RDA was: Cation Exchange Capacity (CEC), Boron (mg/kg), Exchangeable Calcium, Organic Matter, and Exchangeable Magnesium (meq/100 g).
RDA was performed using the rda function from the scikit-bio (v0.6.3) package. We visualised the results with a scatter plot in the 2-dimensional space defined by the main components of the RDA using matplotlib (v3.8.4).
Raw paired-end FASTQ files were reprocessed in QIIME 228 (qiime2-amplicon-2025.7 version) to generate a PICRUSt2-compatible amplicon sequence variant (ASV) table together with the corresponding representative-sequence FASTA. Denoising, was performed with q2-dada2 (v2024.10.0) using denoise-paired. The resulting amplicon sequence variant (ASV) feature table and representative sequences were used to predict functional profiles with PICRUSt229 (v2.4.1). In particular, we predicted Kyoto Encyclopedia of Genes and Genomes (KEGG).
Orthology (KO) terms, Enzyme Commission (EC) numbers, and MetaCyc metabolic pathways for each sample. These functional profiles were aggregated by land-use type; for each functional profile, we generated stacked bar plots showing the relative abundance of the 15 most prevalent functions across land uses. In addition, centered log-ratio (CLR) transformation was applied to each table (after compositional normalization), and principal component analysis (PCA) was performed to visualise patterns of functional differentiation among the four land-use categories. These exploratory analyses were conducted separately for KO, EC, and pathway profiles to compare the degree and structure of functional variation associated with land use. We note that the ASV table resulting from this reprocessing step may differ from the feature table produced by the Illumina pipeline. This approach was adopted because representative sequences were not available from the Illumina pipeline outputs, preventing direct use of that table in PICRUSt2. We note that functional profiles obtained with PICRUSt2 are inferred by reference-genome mapping and phylogenetic placement rather than by direct measurement of functional genes; for this reason, these results should be interpreted cautiously as predictive, also in light of the fact that PICRUSt2 used the QIIME 2-derived ASV table, which may differ from the Illumina pipeline feature table used for other statistical analyses.
