Abstract

Background

   Porous environments constitute ubiquitous microbial habitats across
   natural, engineered, and medical settings, offering extensive internal
   surfaces for biofilm development. While the physical structure of the
   porous environment is known to shape the spatial organization of
   biofilm inhabitants and their interspecific interactions, its influence
   on biofilm community structure and functional diversity remains largely
   unknown. This study employed microfluidic chips with varying
   micropillar diameters to create distinct pore spaces that impose
   different levels of spatial constraints on biofilm development. The
   impact of pore spaces on biofilm architecture, community assembly, and
   metabolic functions was investigated through in situ visualization and
   multi-omics technologies.

Results

   Larger pore sizes were found to increase biofilm thickness and
   roughness while decreasing biofilm coverage over pore spaces. An
   increase in pore size resulted in reduced biofilm community evenness
   and increased phylogenetic diversity. Remarkably, biofilms in 300-μm
   pore spaces displayed the highest richness and the most complex and
   interconnected co-occurrence network pattern. The neutral model
   analysis demonstrated that biofilm assembly within different pore
   spaces was predominantly governed by stochastic processes, while
   deterministic processes became more influential as pore space
   increased. Exometabolomic analyses of effluents from the microfluidic
   chips further elucidated a significant correlation between the
   exometabolite profiles and biofilm community structure. The increased
   community richness in the 300-μm pore space was associated with the
   significantly higher exometabolome diversity.

Conclusions

   Collectively, our results indicate that increased pore space, which
   alleviated spatial constraints on biofilm development, resulted in the
   formation of thicker biofilms with enhanced phylogenetic and functional
   diversity.
   Your browser is not supporting the HTML5 <video> element.
   You may download that video as a file and play it with a player of you
   choice.

   [40]Download video stream.
   [41]Download video file^ (117.2MB, mp4)
   [42]Open in a new tab

   Video Abstract

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s40168-025-02075-0.

   Keywords: Biofilm, Pore space, Spatial constraints, Community assembly,
   Microfluidics, Multi-omics

Background

   Biofilm represents the predominant microbial lifestyle in natural,
   engineered, and medical environments [[43]1]. Porous media, such as
   soil, aquifers, industrial filters, and medical stents, serve as idea
   habitats for biofilm colonization and inhabitation [[44]2–[45]4]. The
   sessile lifestyle confers multifaceted advantages to biofilm
   inhabitants, allowing them to sequester nutrients, withstand shear
   stress, and protect against antimicrobial agents [[46]5]. The biofilm
   development, in turn, exerts profound influences on the physical and
   biochemical environment of porous media systems, playing a crucial role
   in the biogeochemical cycle, bioremediation, and chronic infection
   [[47]6]. Therefore, elucidating the assembly of biofilm communities in
   porous environments is crucial to unveil the pivotal roles of biofilm
   across various ecosystems and develop novel strategies for their
   prevention, management, and engineering.

   The physical structure of the porous environment has been identified as
   one key determinant in shaping microbial communities. The intricate
   physical architecture reduces the transport of nutrients, signal
   molecules, and wastes, generating microscale resource gradients to
   allow a wide microbial diversity via niche specialization [[48]7]. The
   porous environments also enable the coexistence of weak and strong
   competitors by reducing competitive pressure through spatial
   segregation [[49]8–[50]10]. Niche differentiation was thus posited to
   drive the elevated diversity observed in small soil aggregates and
   fine-textured sediments [[51]11, [52]12]. Meanwhile, biofilm formation
   serves as a competitive strategy for space and nutrient [[53]13].
   Previous studies have found that biofilm-forming species eradicate
   biofilm-deficient competitors via smothering them and cutting off
   nutrient access [[54]14, [55]15]. In porous media with limited space
   and restricted nutrient accessibility, biofilm protrudes from the
   grains into the pore space to sequester nutrients and impede the
   colonization of other species on the grains [[56]16]. Consequently, the
   intensified competition driven by biofilm formation contradicts the
   observed enhanced diversity in confined pore spaces, and the underlying
   mechanisms governing the assembly of intricate communities within
   porous media have yet to be fully elucidated. Moreover, the current
   understanding of how pore spaces influence biofilm community assembly
   largely depends on environmental sampling and grain size
   classification, which inherently confound the effects of pore spaces
   with other environmental variables, such as the variations in organic
   and mineral composition [[57]17]. Therefore, the exclusive impact of
   pore spaces on biofilm assembly remains largely unexplored.

   In order to address this knowledge gap, we investigated the assembly of
   biofilm communities and their functional features in pore spaces with
   different levels of spatial constraints (Fig. [58]1). The microfluidic
   chips incorporating arrays of micropillars with varying diameters were
   developed to simulate different porous environments. This approach
   allowed for the in situ monitoring of biofilm structure and morphology
   within the controlled pore spaces. The biofilm community compositions
   were revealed via 16S amplicon sequencing, and the assembly mechanisms
   were elucidated using the Sloan’s neutral community assembly model.
   Exometabolomic analyses were further employed to provide insights into
   the functional properties of biofilms in different pore spaces. Our
   study highlights how expanded pore spaces drive the assembly of
   phylogenetically distant biofilm communities.

Fig. 1.

   [59]Fig. 1
   [60]Open in a new tab

   Microfluidic platform to study soil biofilm assembly and function under
   varying spatial constraints. The microfluidic chips contain
   micropillars of different diameters to impose varying levels of spatial
   constraints on biofilm development. A soil microbial community
   extracted from woodland topsoil was cultivated in the microfluidic
   devices using a soil extract medium prepared from the same soil sample.
   Biofilm functional characteristics were analyzed through exometabolomic
   profiling of the effluent from the microfluidic chips

Results

Reduced spatial constraints promoted biofilm formation

   To investigate the influence of pore space on biofilm assembly, we
   fabricated the microfluidic chips with micropillar arrays featuring
   distinct diameters and spacings at dimensions relevant to biofilm.
   Given that biofilms were generally tens of micrometers thick, pore
   sizes of 20, 50, 100, and 300 μm were chosen to impose different levels
   of spatial constraints (Supplementary Fig. S1). The longitudinal
   section areas were kept the same across different chips, and the height
   of flow channels matched the corresponding micropillar diameter. This
   configuration ensures a uniform porosity of 0.77 for all chips. To
   initiate biofilm development in the porous media, soil microbial
   communities were inoculated into the microfluidic devices and cultured
   using the intensive soil extract medium (ISEM) prepared from the same
   soil sample (Supplementary Fig. S2). Due to the varying volumes of the
   microfluidic chambers with different pore sizes, the fresh medium was
   introduced at a flow rate proportional to each chamber’s volume. This
   approach ensured that the fresh medium remained in the flow chambers
   for an identical duration, subjecting the biofilms to the same flow
   velocity and shear force.

   Biofilms in different pore spaces originated from the small colonies on
   the micropillar surface and displayed distinct morphologies upon
   maturation (Fig. [61]2a & Supplementary Fig. S3). Microscopic image
   analysis was performed to quantify key biofilm characteristics:
   thickness (radial distance from biofilm edges to the micropillar
   surface), roughness (standard deviation of thickness on individual
   micropillars), and coverage (proportion of pore space occupied by
   biofilms) (Fig. [62]2b). Greater biofilm roughness and thickness were
   observed in larger pore spaces with reduced spatial constraints
   (Fig. [63]2c). In the pore sizes ranging from 20 to 100 μm, the
   biofilms extruded from the micropillar surface into the fluid to form
   streamer-like structures, which contributed to a higher biofilm
   coverage compared to that at the pore size of 300 μm (Fig. [64]2a).
   Biofilms in distinct pore spaces also displayed varying compositions of
   extracellular polymeric substances (EPS), with larger pore sizes
   exhibiting a higher proportion of polysaccharides and proteins and a
   decreased proportion of extracellular DNA (eDNA) (Supplementary Figs.
   S4, S5, & S6). This suggests that the larger pore space not only
   enhances biofilm formation by increasing its thickness and roughness
   but also alleviates spatial constraints, leading to a reduced biofilm
   coverage within the pore spaces.

Fig. 2.

   [65]Fig. 2
   [66]Open in a new tab

   Biofilm development in the microfluidic chip. Biofilm morphologies in
   different pore spaces (a). Biofilm cells were stained with DAPI.
   Schematic diagram illustrating biofilm morphology and coverage analysis
   (b). Biofilm thickness was calculated as the average radial distance
   from the biofilm perimeter to the micropillar surface, while roughness
   was quantified as the standard deviation of biofilm thickness on
   individual micropillar surfaces. Coverage represents the proportion of
   pore space occupied by biofilm. The biofilm thickness, roughness, and
   coverage in different pore environments (c). Data are presented as
   mean ± standard deviation. Biofilm morphology and coverage were
   assessed using images sampled from random areas within three
   independent biological replicates per pore space. Twenty-five
   observations per group were analyzed for biofilm thickness and
   roughness, while 15 observations per group were collected for biofilm
   coverage analysis. Different letters indicate significant differences
   (p < 0.05, one-way ANOVA)

Reduced community evenness and increased phylogenetic diversity within larger
pore spaces

   The biofilm communities in different porous media were profiled using
   16S amplicon sequencing. The Shannon diversity index, which considers
   both the number of distinct amplicon sequence variants (ASVs) in the
   community (richness) and their relative abundances (evenness), was
   assessed to evaluate community diversity. The results showed that
   biofilms grown in porous environments exhibited significantly higher
   Shannon diversity and richness compared to those developed in flat flow
   chambers without micropillars (Supplementary Fig. S7). This suggests
   that the porous structure is critical in providing sufficient spatial
   niches to support diverse microbial communities. The Shannon diversity
   of biofilms in porous environments initially decreased as pore size
   increased from 20 to 100 μm but subsequently increased as pore size
   further expanded to 300 μm (Fig. [67]3a). This inverted trend was
   attributed to the reduced evenness and increased richness in larger
   pore spaces, which are positively and negatively correlated with the
   biofilm coverage, respectively (Fig. [68]3b & c). The community
   phylogenetic structure was assessed by calculating the
   abundance-weighted standardized effect size measures of the mean
   nearest taxon distance (weighted SES.MNTD), providing estimates for the
   phylogenetic relatedness within closely-related relatives. The weighted
   SES.MNTD values exhibited a negative correlation with biofilm coverage
   (Fig. [69]3d). Since high biofilm coverage corresponds to high spatial
   constraints in small pore spaces, it indicates that biofilm inhabitants
   under increased spatial constraints are phylogenetically closer and
   exhibit a more uniform distribution.

Fig. 3.

   [70]Fig. 3
   [71]Open in a new tab

   Taxonomic and phylogenetic diversities of biofilms across different
   pore spaces. Shannon diversity of biofilm communities at varying pore
   sizes (a). Negative and positive correlations between community
   richness (b) and evenness (c) with biofilm coverage. Decreased
   abundance-weighted mean nearest taxon distance (weighted SES.MNTD)
   associated with higher biofilm coverage (d). r represents the Spearman
   correlation coefficient, and the shaded area indicates the 95%
   confidence interval. Different letters indicate significant differences
   (p < 0.05, one-way ANOVA, n = 7 independent biological replicates per
   pore size). Relative abundance of whole biofilm communities and
   subcommunities (abundant, intermediate, and rare) at the genus level
   (e)

   The analysis of community composition revealed that Pseudomonas,
   Rhodococcus, and Burkholderia-Caballeronia-Paraburkholderia were the
   most abundant genera across different pore spaces, constituting more
   than 80% of the whole biofilm communities (Fig. [72]3e). To further
   elucidate the difference in biofilms at different pore sizes, we
   compared the abundant (> 0.1%), rare (< 0.01%), and intermediate
   (0.1–0.01%) subcommunities classified based on the relative abundance
   of ASVs in respective porous environments. The Bray–Curtis similarity
   was calculated to quantify microbial composition similarities among
   these subcommunities across varying pore sizes. A total of 42.9% of
   abundant taxa were present across all pore spaces, contributing to a
   Bray–Curtis similarity of 0.7 ± 0.1 (Supplementary Figs. S8 & S9).
   Meanwhile, 85.4% of rare taxa were specific to their respective pore
   space, resulting in a markedly lower community similarity
   (Supplementary Figs. S8 & S9). The increased richness at the pore size
   of 300 μm was attributed to the significantly higher richness of all
   subcommunities compared to those in other pore spaces (Supplementary
   Fig. S10). In comparison to other pore spaces, the abundant
   subcommunities in 300-μm pore space contained more Pseudomonas and less
   Rhodococcus, while the intermediate subcommunities demonstrated higher
   proportions of Bacteroides, Subdoligranulum, and Faecalibacterium
   (Fig. [73]3e).

Stochastic processes drove biofilm community assembly

   We employed the Sloan’s neutral community assembly model to further
   elucidate the biofilm community assembly mechanism in different pore
   spaces. Microbial taxa occurring more or less frequently than expected
   were designated as above or below prediction, respectively. Based on
   the Akaike information criterion, the neutral community model
   outperformed the binomial distribution models, signifying the dispersal
   and ecological drift have a more significant impact than mere random
   sampling of the source metacommunity (Supplementary Fig. S11). The
   biofilm communities at the pore size of 20 μm exhibited the highest
   coefficient of the neutral fit. Increasing pore size led to a decline
   in the neutral model’s explanatory power over community variance (79.9%
   to 31.6%) and the predicted proportion of microbial taxa (79.5 to
   61.4%) (Fig. [74]4a, b, c, d), implying a reduced role of stochastic
   processes in biofilm assembly within larger pore spaces. At the pore
   size of 300 μm, the number of intermediate taxa that occurred more
   frequently than the prediction was increased, whereas more abundant
   ASVs, predominantly assigned to the genera
   Burkholderia-Caballeronia-Paraburkholderia and Pseudomonas, were less
   frequent than expected compared to other pore spaces (Fig. [75]4e, f,
   g, h, i, j, k, l). It suggests that these intermediate taxa possess
   higher migration ability to disperse more widely within the large pore
   space, while more abundant taxa exhibit lower dispersal ability and are
   selected against in the large pore space compared to non-abundant taxa.
   This trend is also corroborated by habitat specialization analysis,
   which reveals a decreased relative abundance of generalists (those
   adaptable to diverse habitats) and a concurrent increase in the
   proportion of opportunists (those falling between generalists and
   specialists) within the larger pore spaces (Supplementary Fig. S12).

Fig. 4.

   [76]Fig. 4
   [77]Open in a new tab

   Fit of the Sloan neutral community model. Each point represents a
   bacterial ASV at the pore size of 20 µm (a, e), 50 µm (b, f), 100 µm
   (c, g), and 300 µm (d, h). Dashed lines illustrate 95% confidence
   intervals around the model prediction. ASVs that occur more frequently
   and less frequently than the prediction are shown in blue and yellow,
   respectively (a, b, c, d). The R^2 demonstrates the goodness of fit to
   the neutral model. The inset pie chart reveals the proportion of ASVs
   below, within, and above neutral expectations. The abundant,
   intermediate, and rare ASVs were colored as orange, green, and purple,
   respectively (e, f, g, h). Relative abundance of different
   subcommunities that fall below, within, and above neutral model
   expectation (i, j, k, l). The numbers in the column denote the count of
   ASVs

   We further investigated the co-occurrence patterns to elucidate the
   potential ecological interactions among biofilm communities in distinct
   pore spaces. The network of taxa at the pore size of 300 μm
   demonstrated the greatest complexity and connectivity because of the
   highest number of nodes and edges (Fig. [78]5, Supplementary Table S1).
   The network topological features, including average path length and
   network diameter, demonstrated reduced values for communities within
   the larger pore sizes (100 μm and 300 μm) compared to those within the
   smaller pore sizes (20 μm and 50 μm) (Supplementary Table S1).
   Additionally, the rare taxa in the pore size of 300 μm exhibited more
   associations between abundant and intermediate subcommunities
   (Fig. [79]5h), indicating a higher level of interconnectedness among
   these subcommunities. Conversely, the network for biofilm communities
   in the pore size of 20 μm exhibited the least nodes and edges, along
   with the lowest clustering coefficient, which quantifies the likelihood
   of two nodes sharing a common neighbor being connected (Supplementary
   Table S1). This suggests that the larger pore space with reduced
   spatial constraints facilitates a wide range of interactions, while the
   biofilm inhabitants within the confined pore space exhibit a low degree
   of association.

Fig. 5.

   [80]Fig. 5
   [81]Open in a new tab

   Network associations of bacterial ASVs in different pore spaces.
   Co-occurrence networks were constructed using Spearman’s rank
   correlation to identify strong and significant correlations (Spearman’s
   r > 0.5 and p < 0.05) (a, b, c, d). The node color represents the
   abundance categories of the ASV, and the node size is proportional to
   its degree. The connection thickness reflects the strength of
   Spearman’s correlation coefficients. A red connection indicates a
   positive correlation, while a green connection denotes a negative
   correlation. The bar graph presents the total number of positive and
   negative connections within each porous environment. Summaries of
   node-connection statistics (e, f, g, h). The black numbers indicate the
   total number of connections within and between subcommunities, and the
   red numbers in parentheses represent the count and percentage of
   positive connections

Increased pore space enhanced function redundancy and exometabolite diversity

   To reveal the functional traits of biofilms assembled in different pore
   spaces, microbial functional potentials were predicted based on the 16S
   amplicon sequencing using Tax4Fun2. At the level 2 of KEGG category,
   genes associated with basic nutritional metabolism and genetic
   information processing, including carbohydrate, nucleotide and glycan
   metabolism, membrane transport, translation, and replication, were
   found to be more abundant in the biofilm at the 20-μm pore size
   compared to those in the larger pore space (Supplementary Fig. S13).
   Meanwhile, the biofilm in 300 μm exhibited the highest relative
   abundance in pathways related to environmental information processing
   and cellular processes, including biofilm development, signal
   transduction, and cell motility (Supplementary Fig. S13). The analysis
   of functional redundancy (FR), based on the average phylogenetic
   distance of community members possessing the same function, revealed a
   larger proportion of functions exhibiting higher FR within the biofilm
   at the 300-μm pore size than those in smaller pore spaces
   (Supplementary Fig. S14). Consequently, the phylogenetically diverse
   biofilm community at the pore space of 300 μm exhibits higher
   multifunctional redundancy.

   As the functional predictions derived from amplicon sequencing do not
   accurately reflect the functions actively expressed, we profiled the
   exometabolites within the effluent from the microfluidic chips via
   metabolomics to elucidate the metabolic functionalities. Procrustes
   analysis revealed a strong correlation between these distinct
   exometabolic profiles and biofilm community structure (Fig. [82]6a).
   The Shannon diversity of the metabolome at the 300-μm pore space was
   significantly higher than those in smaller pore spaces, aligning with
   the increased biofilm community richness (Figs. [83]3b & [84]6b).
   Enrichment analysis was performed to identify metabolic pathways
   significantly enriched within different pore spaces with respect to
   fresh ISEM (Fig. [85]6c). Biofilms within 20-µm pore space exhibited
   the least number of enriched pathways, with a lower enrichment ratio
   compared to those observed in larger pore spaces. Notably, biofilms
   residing in 300-µm pore space displayed unique and significant
   enrichment in pathways associated with D-arginine and D-ornithine,
   butanoate, and alanine, aspartate, and glutamate metabolism, suggesting
   distinct metabolic features for biofilms residing in the largest pore
   spaces (Fig. [86]6c).

Fig. 6.

   [87]Fig. 6
   [88]Open in a new tab

   Enhanced exometabolite diversity in large pore spaces. Procrustes
   analysis revealing the significant correlation between microbial
   communities and exometabolite profiles (a). The Shannon diversity of
   the exometabolites in different pore spaces (b). Enrichment analysis of
   exometabolomic data with respect to fresh ISEM (c). The bubble size
   corresponds to the pathway enrichment ratio for differentially abundant
   metabolites across various pore spaces. Asterisks denote statistically
   significant differences based on the hypergeometric test (***p ≤ 0.001,
   **p ≤ 0.01, *p ≤ 0.05, n = 7 independent biological replicates per pore
   size)

Discussion

   Porous environments serve as the predominant habitats for microbes
   [[89]1, [90]18]. The intricate pore structures within these
   environments modulate diverse microbial interactions and foster the
   development of complex microbial communities [[91]19]. However, the
   role of pore space in shaping biofilm community structure and function
   remains largely unexplored. In this study, we employed
   laboratory-constructed microfluidic chips to investigate the assembly
   mechanisms of soil microbial biofilm communities in pore spaces under
   different levels of spatial constraints. The biofilm in smaller pore
   spaces exhibited reduced thickness and roughness while also
   demonstrating greater coverage by forming a streamer structure within
   these spatially confined pore spaces. Increasing the pore size led to a
   decrease in biofilm community evenness and a concomitant increase in
   phylogenetic diversity. The biofilm inhabiting the 300-μm pore space
   exhibited significantly higher richness compared to those in smaller
   pore spaces. The neutral model analysis revealed a diminishing
   contribution of stochastic processes in biofilm community assembly with
   increasing pore space. Exometabolomic analyses further revealed a
   significant correlation between exometabolite pools and biofilm
   community structure. A significantly higher exometabolome diversity was
   observed in the 300-μm pore space which associated with the increased
   biofilm community richness. Together, these results reveal that large
   pore spaces with reduced spatial constraints promote biofilm formation
   and foster the assembly of a more diverse microbial community, both in
   terms of phylogenetic composition and functional attributes.

   Soil biofilms typically range from 10 to 100 μm in thickness [[92]1,
   [93]20]. In this study, we fabricated microfluidic chips with pore
   sizes ranging from 20 to 300 μm to impose varying degrees of spatial
   constraints on soil biofilm development. According to soil pore size
   classification, pores between 10 and 50 μm were designated as small
   coarse pores, while those higher than 50 μm were defined as wide coarse
   pores [[94]21]. These pore size ranges correspond to mechanical
   fissures, cracks, and biopores (from roots and earthworms), which are
   known hotspots for soil biofilm formation. Previous studies have
   reported divergent effects of pore space on microbial diversity [[95]5,
   [96]22–[97]24]. These varying results were generally attributed to
   complex organic pools, different microbial communities, and distinct
   mineral compositions. To isolate the effect of pore space geometry, we
   investigated the biofilm assembly across different pore spaces with
   identical pore water chemistry, solid substratum, and inoculum
   community composition. Our findings revealed that the smaller pore size
   led to reduced community richness and phylogenetic diversity.
   Consistent with our findings, a previous study in bermudagrass
   ecosystems identified soil texture as the second most influential
   factor (after pH) shaping soil microbial communities, with clay content
   exhibiting negative and positive correlations with bacterial richness
   and evenness, respectively, particularly under wet conditions [[98]24].
   Our analyses further revealed that a stronger fit (R^2) for the Sloan’s
   neutral community model within the smaller pore spaces indicates the
   colonization is predominantly governed by a stochastic, birth–death
   immigration process (Fig. [99]4a, b, c, d) [[100]25]. Deterministic
   processes, such as microbial interactions, appear to exert less
   influence on biofilm assembly in confined porous environments. The
   diminished importance of microbial interaction in shaping biofilm
   community can be attributed to biofilm clogging, which hindered
   nutrient diffusion and metabolite exchange among biofilm inhabitants.
   Consequently, the enhanced microbial diversity within confined pore
   spaces observed in previous studies may arise from factors beyond pore
   size.

   Meanwhile, a significant increase in community richness was observed at
   the 300-µm pore size. Neutral model analysis revealed a higher number
   of abundant taxa were selected against, while more intermediate taxa
   were actively selected at the pore size of 300 μm comparing to smaller
   pore spaces (Fig. [101]4i, j, k, l). It indicates that abundant taxa
   were incapable of occupying all niches within the large pore space,
   thereby enabling the growth of intermediate taxa. This observation
   aligns with the reduced proportion of generalists and the increased
   proportion of opportunists as pore size expands (Supplementary Fig.
   S12). Metabolomic analysis of the 300-μm pore space unveiled greater
   chemical diversity within the exometabolome and the emergence of unique
   metabolic pathways not observed in other small pore spaces, indicating
   extensive metabolic interactions in the large pore space. These
   interactions were further supported by the highest level of network
   complexity and connectivity observed in co-occurrence analysis
   (Fig. [102]5, Supplementary Table S1). It suggests that large pore
   spaces can provide more spatial niches to facilitate the colonization
   of non-abundant taxa and foster the coexistence of phylogenetically
   distant taxa.

   These findings have significant implications for the management,
   control, and engineering of biofilms in porous environments. In
   confined pore spaces, increased biofilm coverage reduces pore
   connectivity, inhibiting biofilm activity and leading to simpler, more
   homogeneous communities. Consequently, biofilms in such environments
   are more susceptible to removal by biocidal agents. In contrast,
   biofilms in larger pore spaces exhibit elevated EPS content, making EPS
   targeting crucial for effective biofilm eradication. Furthermore,
   sufficient pore space supports the development of thicker and more
   metabolically diverse biofilms, suggesting that increasing pore space
   may be an effective strategy for enhancing biofilm functionalities.
   Previous studies have highlighted the superior performance of
   large-pore membrane bioreactors and macroporous biofilm carriers in
   contaminant removal through biocatalysis [[103]26, [104]27].

   This study has elucidated the influence of pore sizes on biofilm
   community assembly and function traits using a combined microfluidic
   and multi-omics approach. Future metagenomic studies could further
   characterize the genomic functional potential of biofilms formed in
   varying pore sizes and at different growth stages. These genomic data
   could corroborate the observed relationship between pore size and
   genomic functional potential, particularly the increased functional
   redundancy within larger pores that may enhance community resilience to
   environmental perturbations, such as fluctuations in pore water and
   nutrient availability, antimicrobial agent exposure, and microbial
   invasion. Further investigation is warranted into the spatial
   organization of different biofilm inhabitants within pore spaces and
   their mutual interactions. The spatial structure of complex biofilm
   communities can be visualized using techniques such as fluorescent gene
   tagging and high-phylogenetic-resolution microbiome mapping by
   fluorescence in situ hybridization (HiPR-FISH) [[105]28]. Additionally,
   the effects of other key physicochemical factors on biofilm development
   within porous ecosystems can be investigated by leveraging microfluidic
   platforms. For example, substratum properties can be assessed by
   coating the internal surfaces of microfluidic devices with polymers or
   minerals through direct adsorption or light curing. Flow conditions can
   be precisely controlled to study hydrodynamic influences, including the
   impact of varying shear forces or unsaturated microenvironments.

Conclusions

   Our study demonstrates that the physical structure of porous
   environments plays a critical role in shaping biofilm community
   structure, assembly processes, and functional diversity. By employing
   microfluidic chips with controlled pore sizes, we observed that larger
   pore spaces increased biofilm thickness and roughness, while decreasing
   biofilm coverage and alleviating spatial constraints. Larger pore size
   resulted in reduced biofilm community evenness and increased
   phylogenetic diversity, with deterministic processes exerting a greater
   influence on community assembly. Exometabolomic analyses of effluents
   from the microfluidic chips further elucidated a significant
   correlation between the exometabolite profiles and biofilm community
   structure. The increased community richness in large pore spaces was
   associated with the significantly higher exometabolome diversity. These
   findings highlight that mitigating spatial constraints in porous
   environments promotes the development of biofilms with greater
   phylogenetic and functional diversity, providing insights into managing
   microbial communities across environmental, industrial, and medical
   contexts.

Methods

Soil sampling

   Soil samples were collected from the woodland topsoil (0–20 cm) at
   Qiyang Red Soil Experimental Station (Hunan province, China, 111°53′E,
   23°27′N) in Mar 2021. The samples were cleaned of stones and roots and
   then homogenized by sieving through a 2-mm mesh.

Soil bacteria extraction and growth medium preparation

   Soil microorganisms were extracted from the collected soil samples
   using Nycodenz gradient centrifugation [[106]29]. Intensive soil
   extract medium (ISEM) was prepared as described in a previous study to
   cultivate soil microbial communities [[107]30, [108]31]. Briefly, 500 g
   of dry soil was subjected to two rounds of extraction with 1.3 L of 80%
   methanol to yield the supernatant, which was subsequently filtered and
   lyophilized. The lyophilized soil extract was redissolved in 200-mL
   water and filtered through a 0.22-μm filter to obtain the soil extract
   solution. ISEM contained 20% soil extract solution, 0.1% mineral salt
   solution, 0.1% vitamin stock solution, 0.2% selenite-tungstate
   solution, and 0.2% trace element solution SL-10. The pH of the ISEM was
   adjusted to 6.8. The ingredient details are available in the
   supplementary information (Supplementary Table S2).

Microfluidic chip fabrication

   Microfluidic chips with arrays of micropillars were fabricated to
   simulate the porous architecture of soil. Four different microfluidic
   chips were designed with micropillars of 20, 50, 100, and 300 µm in
   diameter and height, respectively (Supplementary Fig. S1). The
   corresponding volumes of these microfluidic chambers were 1.2 µL, 3.0
   µL, 6.0 µL, and 18.0 µL. To investigate the influence of the porous
   environment on biofilm development, control flow chambers with
   identical dimensions but without micropillars were also fabricated
   (Supplementary Fig. S15). The direct-write photolithography
   (MicroWriter ML3) and deep reactive ion etching (DRIE) were employed to
   create the microchannel features on silicon wafers. PDMS prepolymer was
   poured over the silicon wafers, degassed under vacuum, and cured in an
   oven at 80 °C for at least 2 h. The structured PDMS was then peeled off
   and bonded to a glass slide after plasma treatment to construct the
   microfluidic devices. The microfluidic chips were sterilized with 75%
   ethanol prior to inoculation.

Biofilm development in microfluidic chips

   The microfluidic chip consisted of a flow channel with two inlets and
   one outlet (Supplementary Fig. S1). One of the inlets was connected to
   a syringe pump to supply fresh ISEM for microbial growth, while the
   other inlet was used for inoculation. The collected soil microbial
   cells were precultivated in ISEM and then resuspended in
   phosphate-buffered saline (PBS) to an OD600 of 0.1. Subsequently, 50,
   125, 250, and 750 µL of soil bacterial suspension were introduced into
   microfluidic chips with micropillars of 20, 50, 100, and 300 µm in
   diameter, respectively. The inoculation port was sealed with PDMS after
   inoculation. Microbial cells were allowed to stand for 2 h to
   facilitate initial attachment to the substratum [[109]30]. Afterward,
   sterile ISEM was introduced into the microfluidic chambers for biofilm
   cultivation at 25 °C in the dark. To maintain a constant retention time
   across different microfluidic chambers, ISEM was then supplied to
   microfluidic chips with micropillars of 20, 50, 100, and 300 µm at flow
   rates of 0.2, 0.5, 1, and 3 μL/min, respectively. Before sampling, the
   tubing and needles at the outlet were replaced with sterile ones, and
   two sterile syringes were connected to the inlet and outlet of
   microfluidic chamber. One milliliter of PBS buffer was repeatedly
   flushed between these two syringes to extract microbial cells in the
   microfluidic chamber. Seven independent biological replicates were
   conducted for each microfluidic design.

Confocal microscopy and image processing

   Biofilms in microfluidic chips were visualized using a confocal
   microscope (AIR, Nikon) after fluorescent staining. The biofilm cells,
   protein, polysaccharide, and eDNA were stained with DAPI, SYPRO Orange,
   Concanavalin A (Con A), and
   7-hydroxy-9H-(1,3-dichloro-9,9-dimethylacridin-2-one) (DDAO),
   respectively. Biofilm thickness and roughness were calculated as
   previously described [[110]30]. The proportional composition of each
   EPS component was determined by dividing its fluorescent area by the
   total fluorescent area of the biofilm cells. Biofilm coverage was
   calculated as the ratio of the biofilm area to the pore area.

16S rRNA amplicon sequencing

   We performed 16S rRNA amplicon sequencing to evaluate the biofilm
   communities assembled in different porous environments. Total DNA was
   extracted with the EZNA Soil DNA Kit (Omega). The V3–V4 regions of
   bacterial 16S rRNA gene were amplified using the primer pair 338F
   (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTA CHVGGGT WTCTAT-3′)
   [[111]32]. The thermal conditions for amplification were as follows:
   initial denaturation at 98 °C for 2 min; 25 cycles at 98 °C for 15 s,
   55 °C for 30 s, and 72 °C for 30 s; and final elongation at 72 °C for
   30 s. The PCR amplicons were purified with VAHTS™ DNA Clean Beads
   (Vazyme) and then sequenced on the Illumina MiSeq PE300 platform.

Exometabolomic analysis

   Prior to extracting microbial communities from microfluidic devices,
   100 μL of effluent was collected at the outlet for metabolomic
   analysis. The effluent was passed through a 0.2-μm filter to remove
   bacterial cells. The filtrate was then vortexed with 400 μL of
   pre-cooled methanol/acetonitrile/water solution (2:2:1, v/v/v),
   incubated at − 20 °C for 10 min, and centrifuged at 14,000 g at 4 °C
   for 20 min. The supernatant was dried in vacuum and redissolved in 100
   μL of acetonitrile aqueous solution (acetonitrile:water = 1:1, v/v).
   After centrifuged at 14,000 g at 4 °C for 15 min, the supernatant was
   subjected to LC–MS/MS analysis. Quality control (QC) samples were
   prepared by combining equal aliquots from all effluent samples and
   measured before analysis and after every 10 runs to evaluate the
   analytical variance. The metabolomic features were characterized using
   ultrahigh performance liquid chromatography-tandem time-of-flight mass
   spectrometry (UHPLC-Q-TOF MS). The ESI source temperature was set to
   600 °C, with the ion spray voltage floating at 5500 V in both positive
   and negative modes. TOF MS scan m/z range was 60–1000 Da, and the
   product ion scan m/z range was 25–1000 Da.

Metabolomic data analysis

   ProteoWizard was used to convert the raw data into mzXML format, which
   was then processed by XCMS for peak alignment. XCMS was configured with
   the following settings: method = “centWave”, ppm = 10, peakwidth = c
   (10, 60), prefilter = c (10, 100), mzwid = 0.025, minfrac = 0.5, and
   bw = 5. Only metabolites at level 2 identification (putatively
   annotated compounds) were used for downstream analysis. The
   differentially abundant metabolites (DAMs) in the effluents, relative
   to ISEM, were identified through OPLS-DA and t-test (VIP > 1,
   p < 0.05). Metabolic pathway enrichment analysis was carried out using
   the MetaboAnalyst 5.0 web-based platform ([112]www.metaboanalyst.ca) to
   determine the KEGG metabolic pathways significantly enriched in the
   DAMs [[113]33].

Bioinformatic analysis

   The raw amplicon sequencing data were processed using the DADA2 plugin
   in QIIME2 ([114]https://view.qiime2.org) [[115]34, [116]35]. The
   Cutadapt plugin was employed to trim the primers (cutadapt trim-paired)
   [[117]36]. The sequences were then denoised using DADA2 based on the
   quality score, clustered into amplicon sequence variants (ASVs), and
   rarefied to the minimum number of reads per sample (21,247 reads). The
   ASV sequences were aligned to a pre-trained SILVA database using the
   QIIME2 feature-classifier plugin (release 138, bacteria) [[118]37].
   Richness and Shannon diversity indices were calculated using the vegan
   package [[119]38]. Following the classification in the previous studies
   [[120]39], ASVs were categorized as abundant (> 0.1%), rare (< 0.01%),
   or intermediate (0.01–0.1%) based on their relative abundances in each
   porous media. The functional profiles of microbial communities were
   predicted by Tax4Fun2 [[121]40]. Co-occurrence networks were
   constructed using Spearman’s rank correlation coefficients between ASVs
   present in at least two samples. Only robust correlations (|r|> 0.5)
   with statistically significant (p < 0.05) were included in the
   networks. The networks were visualized using the interactive platform
   Gephi [[122]41–[123]43]. Network topological features, including the
   number of nodes and edges, modularity, average clustering coefficient,
   network diameter, average path length, and average degree, were
   analyzed using Gephi. Phylogenetic diversity was evaluated through the
   abundance-weighted standardized effect size measures of the mean
   nearest taxon distance (weighted SES.MNTD), using the picante package
   [[124]44].

Community assembly analysis

   The Sloan neutral community model (NCM) was employed to elucidate the
   assembly mechanisms of communities in different porous environments
   [[125]25, [126]45]. The relationship between ASV occurrence frequency
   and their relative abundance across the wider metacommunity was
   predicted to determine the contribution of stochastic process. The
   goodness of fit to the neutral model was assessed using R^2 through
   nonlinear least-squares fitting. Bacterial taxa exceeding, aligning
   with, or falling below the model-predicted expectations were identified
   by comparing their frequencies with the 95% confidence interval limits
   of this best-fit neutral model. The fit of the neutral model was
   compared with that of the binomial distribution model to determine the
   effects of ecological drift and dispersal limitation on community
   assembly [[127]46].

Habitat specialization analysis

   The assessment of habitat specialization was conducted using Levins’
   niche breadth [[128]47]. Levins’ niche breadth was calculated as
   follows:
   [MATH:
   <mrow><msub><mi>B</mi><mi>j</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow
   ><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</m
   i></msubsup><msubsup><mi>P</mi><mrow><mi
   mathvariant="italic">ij</mi></mrow><mn>2</mn></msubsup></mrow></mfrac><
   /mrow> :MATH]

   where
   [MATH: <msub><mi>B</mi><mi>j</mi></msub> :MATH]
   indicates the niche breadth of ASV
   [MATH: <mi>j</mi> :MATH]
   , N is the total number of communities, and
   [MATH: <msub><mi>P</mi><mrow><mi
   mathvariant="italic">ij</mi></mrow></msub> :MATH]
   is the proportion of ASV
   [MATH: <mi>j</mi> :MATH]
   in community i. Niche breadth was calculated based on the relative
   abundance of ASVs within each porous environment. The occurrences of
   ASVs wer generated by the simulation of 1000 permutations using the
   quasiswap algorithm of the EcolUtils R package [[129]48]. ASVs are
   defined as generalists if their observed occurrence surpasses the upper
   limit of the 95% confidence interval. Conversely, any ASV with observed
   occurrence below the lower 95% confidence interval is categorized as a
   specialist. ASVs with occurrences within the 95% confidence interval
   are considered opportunist taxa [[130]49].

Sample size determination

   Previous research has suggested using more than five biological
   replicates in amplicon sequencing to detect differences in microbial
   communities effectively [[131]50]. In this study, we used seven
   biological replicates in our amplicon sequencing experiments to enhance
   data comparability and achieve higher resolution in differentiating the
   effects of varying pore sizes. To evaluate the significance of pore
   space effects on community structure across different biological
   replicates, we conducted three nonparametric analyses: nonparametric
   multivariate ANOVA using distance matrices (Adonis), analysis of
   similarity (ANOSIM), and multiresponse permutation procedure (MRPP).
   The Bray–Curtis similarity index was used to construct the distance
   matrix for all three analyses. Statistical significance was assessed
   using Monte Carlo permutations. All analyses were conducted using the
   vegan package [[132]38]. By comparing community structures across
   different numbers of biological replicates, we confirmed that seven
   replicates were sufficient to distinguish differences in community
   structures (Supplementary Table S3). Furthermore, we conducted a power
   analysis to determine the appropriate sample size for imaging analysis.
   Using a significance level of 0.05, a statistical power of 0.8, and an
   effect size of 0.5 across four groups, the power analysis indicated
   that 12 observations per group were necessary. Based on these results,
   we collected 15 observations per group for EPS and biofilm coverage
   analysis and 25 observations per group for biofilm thickness and
   roughness analysis from random areas within three biological
   replicates.

Supplementary Information

   [133]40168_2025_2075_MOESM1_ESM.docx^ (25.7MB, docx)

   Supplementary Material 1: Supplementary Fig. S1. Microfluidic chips
   with micropillar arrays designed to simulate porous environments. The
   chip features micropillars with diameters of 20, 50, 100, and 300 µm,
   with corresponding heights. It includes one outlet, one inoculation
   port, and one medium inlet. The inoculation port is positioned
   downstream of the medium inlet to prevent contamination. The pore size
   is defined by the minimum spacing between adjacent micropillars (a).
   Top views (b-e) and side views (f-i) of microfluidic chips with
   micropillars of varying diameters. Microscopic top-down views (j-m) and
   scanning electron microscopy (SEM) images (n-q) of the micropillar
   arrays. Supplementary Fig. S2. Flowchart illustrating the experimental
   design. A soil microbial community was extracted and cultured in
   microfluidic chips with varying pore sizes using a soil extract medium
   derived from the same soil sample. Following biofilm development,
   biofilm morphologies were imaged, microbial community compositions were
   analyzed via 16S amplicon sequencing, and functional properties were
   assessed using exometabolomics. Supplementary Fig. S3. The development
   of biofilm architecture in the porous environment of various sizes: 20
   μm (a), 50 μm (b), 100 μm (c), and 300 μm (d), respectively.
   Supplementary Fig. S4. 3D visualization of biofilms grown in
   microfluidic chips with varying pore sizes. Biofilms were stained with
   DAPI (blue) for cells, SYPRO Orange (green) for proteins, Concanavalin
   A (Con A, red) for polysaccharides, and
   7-hydroxy-9H-(1,3-dichloro-9,9-dimethylacridin-2-one) (DDAO, purple)
   for eDNA. Supplementary Fig. S5. Orthogonal views of biofilms grown in
   the microfluidic chips with different pore sizes. Supplementary Fig.
   S6. The relative abundance of different EPS components in biofilms.
   Extracellular proteins, polysaccharides and eDNA were stained with
   fluorescent dyes. The relative abundance of each EPS component was
   determined as the ratio of its fluorescence area to the total
   fluorescence area of biofilm cells. 15 observations per group were
   collected for EPS analysis. Different letters indicate significant
   differences (p < 0.05, one-way ANOVA). Supplementary Fig. S7. Biofilms
   developed in microfluidic chips with porous structures exhibit
   significantly higher Shannon diversity (a), richness (b), and evenness
   (c) compared to those in flat chambers. Different letters indicate
   significant differences (p < 0.05, one-way ANOVA). Supplementary Fig.
   S8 Pairwise comparison of Bray-Curtis similarity between different
   porous environments for abundant (a), intermediate (b) and rare taxa
   (c). Significance of results was evaluated using ANOVA and indicated by
   different letters. Supplementary Fig. S9 Venn diagram illustrating the
   numbers of unique and shared ASVs among various porous environments for
   different subcommunities: Abundant (abundant taxa), Intermediate
   (intermediate taxa), Rare (rare taxa), and Whole (all bacterial taxa).
   Supplementary Fig. S10 The taxonomic richness diversity across
   different porous environments for abundant, intermediate, and rare
   taxa. Statistical significance was assessed using ANOVA, with different
   letters indicating significant differences. Supplementary Fig. S11
   Comparison of the maximum likelihood fit between the neutral and
   binomial models. Supplementary Fig. S12 Relative abundance of
   generalists, opportunists and specialists in respective porous
   environments. Significant differences were evaluated using ANOVA and
   denoted with different letters. Supplementary Fig. S13 Relative
   abundances of KEGG level2 pathways. Different letters indicate
   significant differences (p < 0.05, one-way ANOVA). Supplementary Fig.
   S14 Differences in functional redundancy in different pore
   environments. 300 µm vs 20 µm (a); 300 µm vs 50 µm (b); 300 µm vs 100
   µm (c). A log ratio greater than zero denotes greater functional
   redundancy in the 300 µm pore environment. Supplementary Fig. S15 The
   photograph of control microfluidic chips without micropillar arrays.
   Top view of the microfluidic chips (a-d). Cross-sectional view of the
   flat chambers with varying heights (e-h). Supplementary Table S1.
   Network-level topological features. Supplementary Table S2. Recipes for
   soil extract solutions. Supplementary Table S3. Significance tests of
   pore space effects on the biofilm community structure across different
   biological replicate numbers and statistical approaches.

Acknowledgements