The soil sampling sites were located at Linyi City (34°22′ ~ 36°13′N, 117°24′ ~ 119°11′E) of Shandong Province, Eastern China (Figure 1A). As a national manufacturing base in China, Linyi City has a variety of industrial production categories. This area experiences a temperate monsoon climate with a mean annual precipitation of 799.9 mm (~2.2% occurs during October to December) and a mean annual temperature of 15.0 °C.

To capture the natural heterogeneity of industrial environments and robustly test the generality of the cement-hardening effect on soil ecosystems, a stratified selection approach was employed. We selected five factory sites representing two major industrial categories (chemical manufacturing and synthetic materials & processing) in the region, each defined by distinct production activities and associated potential soil contamination profiles. This design ensured our findings were not idiosyncratic to a single factory type but reflected broader patterns across industrial landscapes. All sites were located within Linyi City to control for macro-climatic and broad soil-type variations. Critically, within each factory, the paired plot design was used to isolate the effect of cement hardening from pre-existing land-use history. The selected factories were: Fenyan Co., LTD. (FY, 35°28′N, 117°44′E), which produced chemical fertilizer in Pingyi County, representing the chemical manufacturing category; Zhuoyue Wood Industry Co., LTD. (ZY, 35°12′N, 118°23′E, plywood) and Anfen Chemical Co., LTD. (AF, 35°08′N, 118°26′E, polyethylene) in Lanshan and Luozhuang Districts, respectively, representing the synthetic materials and processing category; Yaocheng Nickel-Chromium Alloy Co., LTD. (YC, 35°06′N, 119°09′E, nichrome) in Ju′nan County, also under chemical manufacturing with a focus on metallurgy; and Shenyan Group Co., LTD. (SY, 34°54′N, 118°18′E, zetar) in Luozhuang District, representing synthetic materials & processing.

Soil sampling was conducted in November 2023 across five factories in Linyi City using a paired sampling design to isolate the effect of cement hardening. Within each factory, six pairs of adjacent sampling plots (5 m × 5 m each) were established, uniformly distributed across the factory grounds. Each pair consisted of one cement-hardened (CH) plot and one adjacent unsealed bare soil (CK) plot (within 1 ~ 3 m), both sharing identical pre-sealing land-use history. For CH plots, the cement surface debris was first removed by brushing and rinsing with distilled water to minimize surface contamination of subsurface sample. A crawler-mounted core drill (SE3000, Qianhe Environment, China) cooled with sterile distilled water was then used to cut through the cement layer and extract a 20-cm diameter cement core. After manual removal of the cement core, the underlying undisturbed soil surface was exposed. A sterilized polyethylene tube (6.0 cm diameter) was inserted vertically into the soil profile to collect a subsurface soil core from 0 ~ 20 cm depth. For adjacent CK plots, soil cores (0 ~ 20 cm depth, 6.0 cm diameter) were collected vertically using a manual T-type soil sampler. The sampler was sterilized with 75% ethanol between plots. In each plot, three soil cores were collected from random locations and homogenized to form one composite sample per plot. A total of 60 soil samples were obtained, comprising 30 CK samples and 30 CH samples. All soil samples were passed through a 2 mm mesh to remove stones and plant residues, and were brought back to the laboratory. A portion of the soil samples was air-dried for soil physical–chemical analysis, and the second portion was stored at −80 °C for soil microbial DNA extraction.

Soil moisture was determined gravimetrically by comparing the weights of fresh soils before and after ovendrying at 105 °C until constant weight. Soil pH and electrical conductivity (EC) values were measured using a soil-to-water ratio of 1:2.5 and 1:5 (w/v) suspension by a pH meter (BPH-220, Bell Instrument, China) and a conductometer (DDSJ-308F, Leici Instrument, China). Soil organic matter (OM) was estimated using the loss on ignition method. Samples were dried at 105 °C for 24 h prior to weighing, heated at 550 °C for 6 h and reweighed to determine the loss of OM by weight as a percentage. Total carbon and nitrogen contents were determined by the Vario EL III elemental analyzer (El Vario analyzer, Elementar, Germany). Soil total phosphorus, total potassium, and total heavy metals (chromium, nickel, plumbum, and arsenic) contents were determined using the inductively coupled plasma-optical emission spectrometer (ICP-OES) technique (Optima 5300DV, PerkinElmer, USA) after digesting by an acid mixture of nitric acid (HNO₃), hydrofluoric acid (HF), and perchloric acid (HClO₄) (10:5:2, v/v) using the electric heating plate (XJS20-42, Laboratory Instrument, China). It was noted that the pretreatment procedures differed among analyzes according to the specific property measured: state variables such as moisture, pH, and EC were analyzed on fresh soils; total elemental contents required acid digestion; and organic matter was determined via combustion. All samples for any given analysis were subjected to the identical pretreatment and analytical protocol to guarantee comparability across the CH and CK groups.

Total genome DNA was extracted from 500 mg of frozen soil samples using CTAB-SDS method (Zhou et al., 1996). DNA quality and purity was checked with 260/280 and 260/230 nm wavelength ratios using the NanoDrop™ (Thermo Scientific, Massachusetts, USA). The primer sets of 341F (5’-CCTAYGGGRBGCASCAG-3’) and 806R (5’-GGACTACNNGGGTATCTAAT-3’) were used to amplify the V3 ~ V4 region of bacterial 16S rRNA gene (Berg et al., 2012). The primer pairs of ITS1-1F-F (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS1-1F-R (5’-GCTGCGTTCTTCATCGATGC-3’) were selected to amplify the ITS1-1F region of fungal rRNA gene (Li et al., 2021). The 30 μL polymerase chain reaction (PCR) mix consisted of 13 μL sterilized H₂O, 15 μL High-Fidelity PCR Master Mix (New England Biolabs, USA), 0.5 μL each primer, and 1 μL template DNA (10 ng). DNA samples were amplified using the following thermal cycling conditions in three replicate runs: initial denaturation at 98.°C for 1 min, followed by 30 cycles of denaturation at 98.°C for 10 s, annealing at 50.°C for 30 s, and elongation at 72.°C for 60 s, and a final extension step at 72.°C for 5 min. Duplicate PCR products mixed with the same volume of 1 × loading buffer (contained SYB green) were pooled and quality checked on 2% agarose gel. The amplicons with bright main strip between 400 ~ 600 bp were chosen and purified using the E. Z. N. A.™ Gel Extraction Kit (Omega Bio-Tek Inc., USA). Sequencing libraries were generated and library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific, USA) and Agilent Bioanalyzer 2,100 system (Agilent Technologies, Palo Alto, USA). At last, the library was sequenced on an Ion S5™ XL platform (E-gene, Shenzhen, China) and 400 bp/600 bp single-end reads were generated.

Microbiome bioinformatics were performed with quantitative insights into microbial ecology 2 (QIIME 2) with slight modification according to the official tutorials. Briefly, raw sequence data were demultiplexed using the demux plugin following by primers cutting with cutadapt plugin (Martin, 2011). Sequences were then quality filtered, denoised, merged and chimera removed using the divisive amplicon denoising algorithm 2 (DADA 2) plugin (Callahan et al., 2016). Non-singleton amplicon sequence variants (ASVs) were aligned with mafft (Katoh et al., 2002) and used to construct a phylogeny with fasttree2 (Price et al., 2010). Taxonomy was assigned to ASVs using the classify-sklearn naïve Bayes taxonomy classifier in feature-classifier plugin (Bokulich et al., 2018) against the Greengenes (McDonald et al., 2012) and UNITE (Abarenkov et al., 2010) databases for bacterial and fungal sequences.

Microbial co-occurrence networks were constructed by the Molecular Ecological Network Analyzes pipeline (MENAP, http://ieg2.ou.edu/MENA/) based on the Spearman’s correlation matrices (Wen et al., 2022). The correlation thresholds of r threshold > 0.6 and p threshold < 0.05 were used to build the four networks of soil bacterial and fungal communities from CK and CH samples. Network topological properties, including connectance, average degree, average path length, diameter, mean clustering coefficient, Centralization, and modularity were calculated (Deng et al., 2012). The robustness was performed to evaluate the community stability of CK and CH network. In this study, a certain proportion of nodes was randomly removed to simulate random species removal, and certain numbers of module hubs were removed to simulate targeted removal. The robustness of the network was measured by the proportion of nodes remaining after deletion of the specified species.

Niche breadth was used to measure habitat specialization patterns among ASVs using the Spaa R package. The higher niche breadth values of ASVs represented taxa that were widely distributed on a large scale with relatively even abundance, while those with lower values indicated taxa occurred in fewer habitats with uneven distribution patterns. The niche breadth values of CK and CH soil samples were calculated based on species abundance and evenness.

ASVs occurrences were simulated to determine the niche specialization through 1,000 permutations using quasiswap permutation algorithms implemented in the EcolUtils R package. Generalist species typically exhibited broader fundamental niches compared to specialists. In this study, ASVs were classified as generalists or specialist if their observed occurrence exceeded the upper or fell below the lower 95% confidence interval. ASVs with niche breadth values within the 95% confidence interval range were categorized as neutral taxa.

To assess the importance of deterministic mechanisms from stochastic mechanisms in microbial community assembly, the β nearest taxon index (βNTI) was calculated based on taxonomic metrics using the NST R package. Homogeneous selection (βNTI < −2) and heterogeneous selection (βNTI > + 2) were used to describe the deterministic assembly process of CK and CH soil microorganisms.

Potential bacterial functions among soil microbiota in CK and CH samples were predicted by BugBase (Ward et al., 2017), phylogenetic investigation of communities by reconstruction of unobserved states 2 (PICRUSt2) (Douglas et al., 2020), and FAPROTAX (Louca et al., 2016) softwares, and fungal functions were predicted based on fungal taxa using the FUNGuild database (Nguyen et al., 2016). BugBase provided an algorithm that predicted organism-level coverage of functional pathways as well as biologically interpretable phenotypes. PICRUSt2 contained an updated and larger database of gene families and reference genomes, provided interoperability with any denoising algorithm. PICRUSt2 inferred metagenomic potential based on phylogenetic placement and reference genomes, which might be affected by horizontal gene transfer and database gaps. FAPROTAX was a database that mapped prokaryotic clades (e.g., genera or species) to established metabolic or other ecologically relevant functions, using the current literature on cultured strains. FUNGuild was an open annotation tool for parsing fungal community datasets by ecological guild. FUNGuild assigned ecological roles based on published literature, but fungal nutritional modes were often facultative, and assignments were taxonomy-dependent. Firstly, the rarified taxon tables were first translated into function tables, based on taxon-function annotations in the BugBase, PICRUSt, FAPROTAX, and FUNGuild database. Then the relative abundance of each functional group was calculated based on the number of sequences per sample.

The alpha-diversity indices of Chao1 and Shannon index were calculated to measure bacterial and fungal community structures, and richness, Pielou’s evenness, and Shannon index were used to estimate bacterial and fungal community functions. Non-metric multidimensional scaling (NMDS) ordination and principal component analysis (PCA) based on the Bray−Curtis distance metrics were used to present the heterogeneities of microbial community structures and functions between the CK and CH soil samples. The Anosim analysis was performed to show the difference significance. The contributions of species replacement and species richness difference on the total beta-diversity of microbial community were analyzed using adespatial R package (Shen et al., 2020). The significant differences in soil physicochemical factors, alpha-diversity indices of microbial community structure and function, and relative abundances of dominant microbial phyla, genera, and functional traits between the CK and CH soils were analyzed by Wilcoxon signed-rank test. The significant effects of factory location and cement-hardened time on each soil property were estimated by One-way ANOVA. Canonical correlation analysis (CCA) and Mantel tests were used to link the microbial community structures to soil variables. The soil variables significantly (p < 0.05) associated with microbial community changes were selected for CCA plotting. Pearson correlation analysis was used to test the interrelation coefficients of dominant bacteria and fungi at the phylum and genus levels and soil environmental variables. Partial least square path modeling (PLSPM) using packages amap, shape, diagram, and plspm was constructed to analyze the effects of cement hardening on soil variables, microbial community composition and diversity, and microbial function composition and diversity.

Soil physicochemical properties exhibited considerable variations across the five factory sites. Soil variables from the Yaocheng Nickel-Chromium Alloy Co., LTD. (YC, nichrome production) and Fenyan Co., LTD. (FY, chemical fertilizer production) differed visibly from those of the other three factories (Supplementary Figure S1). Notably, the YC site showed the lowest levels of soil organic matter (OM), total carbon (TC), and total nitrogen (TN), but the highest concentrations of total nickel (TNi) and total chromium (TCr), which reflected the inherent regional variations in industrial history and baseline conditions (Table S1). The soil physicochemical properties under cement hardening (CH) tended to deviate from those of the bare soils (CK) as indicated by principal component analysis (PCA). In detail, cement hardening significantly (p < 0.05) increased soil moisture by 45.2% and decreased soil contents of OM, TC, and total potassium (TK) by 3.8 g/kg, 10.0 g/kg, and 1.3 g/kg, respectively (Table 1). The significant (p < 0.05) effects of factory location and cement-hardened time on soil electrical conductivity (EC) value and soil nutrients such as OM, TC, TN, total phosphorous (TP), and TK were observed. Soil EC and pH values were higher in CH soils than that of CK soils, with an increase of 95.1 μS/cm and 0.2 units. In addition, the impact of cement hardening on soil TN was consistent with TC, but increased soil available nitrogen (AN) level by 16.6%. Although no significant sealed effects were observed on total heavy metal levels, cement hardening decreased the soil contents of TCr, TNi, total plumbum (TPb), and total arsenic (TAs).

There was no significant (p > 0.05) difference in bacterial Chao1 value between CK and CH soils (Figure 2A), but the Shannon indexes decreased as the soils were sealed with cement (Figure 2B). Both the fungal richness and diversity of the CH soils were significantly (p < 0.005) lower than those of the CK soils (Figure 2B). Non-metric multidimensional scaling (NMDS) plots revealed the divergent bacterial and fungal communities after the soil seal (Figures 2C,D), and the Anosim tests further indicated the bacterial (R = 0.627, p < 0.001) and fungal (R = 0.245, p < 0.001) community structures of the CK soils were highly dissimilar to that of the CH soils. The beta-diversity decomposition analysis showed that the species replacement processes dominated the bacterial (74.9%) and fungal (71.1%) community compositional dissimilarities among all study samples, while the richness difference processes only contributed 25.1 and 28.9% on average. Cement hardening significantly (p < 0.05) decreased the relative abundances of Actinobacteria, Proteobacteria, Acidobacteria, Gemmatimonadetes, Bacteroidetes, Verrucomicrobia, and Planctomycetes in the bacterial community, but increased (p < 0.001) the relative abundance of Firmicutes in the CK soils (Figure 2E). For the fungal community, the decreased relative abundances under cement seal were observed for Basidiomycota, Zygomycota, Glomeromycota, and Chytridiomycoa, while increased the relative abundance of Ascomycota. At the genus level, cement hardening led to a significant (p < 0.05) decrease in the relative abundances of bacterial genera Sphingomonas, RB41, MND1, and Nocardioides, as well as fungal genera Fusarium, Mortierella, and Humicola. Conversely, it significantly (p < 0.001) increased the abundances of bacterial genera including Bacillus, Ochrobactrum, Paenisporosarcina, Paenibacillus, and Brevibacillus.

Soil cement hardening distinctly changed the topological properties of soil microbial co-occurrence networks, especially for the fungal network (Table 2). The bacterial network of CH soils was greatly larger (e.g., more nodes, links, and diameter) than that of CK soils. Although the nodes numbers of the fungal network were slightly decreased from 382 (CK) to 375 (CH), the links number and diameter were 2.3 times and 1.2 times higher that of CK soils, respectively. In addition, more complex (e.g., higher connectance, average degree, and average clustering coefficient) CH soil networks of both bacterial and fungal communities were observed compared to CK soils. Conversely, the modularity net and modularity random of the microbial networks were decreased with the cement seal of the soil.

Both random and targeted taxa removal were simulated to evaluate the robustness of microbial networks, quantified as their resistance to node loss (Figure 3). Under random taxa loss, bacterial and fungal networks in cement-sealed soils exhibited robustness comparable to those in bare soils (Figures 3A,C). Notably, based on targeted taxa removal, the bacterial network demonstrated higher robustness in CH soils compared to CK soils (Figure 3B). In contrast, fungal networks showed an inverse response pattern, with CH soils demonstrating markedly reduced robustness relative to CK soils (Figure 3D).

The niche breadth was applied to assess microbial habitat specialization patterns following cement hardening. Cement hardening distinctly reduced the niche breadth of bacterial community (Figure 4A), accompanied by a marked increase in specialist species proportion and a corresponding decrease in generalist species (Figure 4B). In contrast, the fungal community exhibited only marginal elevation in niche breadth under cement-hardened conditions. The percentage of generalist fungal species increased from 2.8% (CK) to 6.5% (CH), but neutral species remained the dominant component of the fungal community in both CK and CH soils. Deterministic processes exerted a stronger influence on the assembly of both bacterial and fungal communities in CH soils compared to CK soils (Figures 4C,D). Cement hardening significantly enhanced heterogeneous selection process in microbial community assembly. Particularly for fungi, the deterministic heterogeneous selection process accounted for nearly all (96%) of the observed phylogenetic variation in CH soils, which presented a larger extent than observed in bacteria.

The putative metabolic or ecological functions of bacterial and fungal amplicon sequence variants (ASVs) were comparatively analyzed between CK and CH soils. Bacterial α-diversity (richness) remained unaffected by cement hardening, while both the Pielou evenness and Shannon diversity indices exhibited significant reductions (p < 0.001) (Figure 5A). In contrast, fungal functional α-diversity demonstrated a pronounced decline across all metrics, including richness, Pielou evenness, and Shannon diversity, in response to cement hardening (Figure 5B). PCA revealed that significant divergences (p < 0.001) in the functional composition of both bacteria and fungi were observed under cement-hardened conditions (Figures 5C,D).

As predicted by the Bugbase analysis, soil cement hardening significantly decreased the relative abundances of bacterial functional groups associated with forms biofilms, Gram⁻ bacteria, and both aerobic and anaerobic taxa, while enriched potentially pathogenic bacteria, mobile elements, Gram⁺ bacteria, and facultatively anaerobes (Figure. 5E). Phylogenetic investigation of communities by reconstruction of unobserved states 2 (PICRUSt2) analysis suggested that cement hardening exerted significant negative effects (p < 0.001) on pathways related to membrane transport, glycan biosynthesis and metabolism, and xenobiotics biodegradation and metabolism. In contrast, functional groups associated with transcription, nucleotide metabolism, cell motility, terpenoid and polyketide metabolism, cofactor and vitamin metabolism, and amino acid metabolism were markedly enhanced. FAPROTAX annotation demonstrated a pronounced decline (p < 0.001) in chemoheterotrophy, aerobic chemoheterotrophy, and aromatic compound degradation in CH soils compared to CK soils. Conversely, metabolic pathways linked to fermentation, nitrate reduction, sulfur compounds/sulfate respiration, and methanogenesis exhibited an increasing trend. Fungal functions were also significantly (p < 0.05) impacted, with cement hardening reducing the relative abundances of dominant ecological roles such as endophyte, wood/soil/litter/plant/dung saprotrophs, animal pathogen, lichen parasite, fungal parasite, ectomycorrhizal fungi, and epiphyte across all sampling sites (Figure 5F).

Mantel test revealed that soil moisture, total carbon, and factory location exerted a significant (p < 0.05) influence on both bacterial and fungal communities (Figure 6A). Additionally, total Pb, total nitrogen, and organic matter content had strong effects on bacterial community composition, whereas total Cr and total Ni significantly (p < 0.05) shaped fungal community structure. Canonical correspondence analysis (CCA) further demonstrated that bacterial and fungal community structures were strongly correlated with soil moisture, electrical conductivity, factory location, and the total levels of Cr, Ni, Pb, nitrogen, and carbon (Figures 6B,D). Furthermore, cement-hardened time, total potassium, and organic matter also significantly (p < 0.05) affected the fungal community structure. A heatmap analysis illustrated the influence of soil environmental variables on the relative abundance of dominant bacterial and fungal phyla (Figures 6C,E). Specifically, soil moisture, electrical conductivity, and pH were negatively correlated with the abundances of Acidobacteriota and Actinobacteriota but positively associated with Nitrospirae. Total heavy metal contents in soil positively affected the relative abundances of Acidobacteria, Actinobacteria, Gemmatimonadetes, and Bacteroidetes. Furthermore, Acidobacteria and Gemmatimonadetes exhibited positive correlations with total nitrogen, total carbon, and organic matter. For the fungal phyla, soil moisture and pH significantly affected the relative abundance of dominant fungi. Notably, heavy metal contents contributed more substantial role than soil nutrient components such as available nitrogen and total potassium in driving fungal distribution differences between bare soil and cement-hardened soil. At the genus level, soil moisture, total carbon, and factory location were significantly (p < 0.05) correlated with the bacterial abundances of Ochrobactrum, Nocardioides, and Meiothermus. Total plumbum and arsenic contents negatively affected the relative abundances of Ochrobactrum, RB41, MND1, and Gaiella. Meiothermus, Bacillus, Brevibacillus, and Paenibacillus exhibited positive correlations with organic matter, total phosphorous, and available nitrogen. For the fungal genus, Mortierella was positively correlated with soil chromium, total nickel, and total potassium, but negatively correlated with soil moisture and electrical conductivity. Furthermore, soil organic matter and total carbon showed positive correlations with Alternaria and Phoma.

Partial least squares path modeling (PLSPM) analysis elucidated the effects of cement hardening on the diversity of bacterial and fungal community structure and function (Figures 7A,B). Cement hardening exerted a direct influence on bacterial community composition and diversity, which in turn indirectly affected bacterial functional composition and diversity. Additionally, soil properties and nutrient levels exhibited significant relationships with both bacterial community diversity and functional diversity. In contrast, cement hardening directly impacted fungal community composition and functional diversity, while indirectly modulated fungal functional composition and community diversity. Furthermore, soil properties and nutrients played a direct role in shaping fungal functional diversity.