Ecosystem size-induced environmental fluctuations affect the temporal dynamics of community assembly mechanisms

Understanding processes that determine community membership and abundance is important for many fields from theoretical community ecology to conservation. However, spatial community studies are often conducted only at a single timepoint despite the known influence of temporal variability on community assembly processes. Here we used a spatiotemporal study to determine how environmental fluctuation differences induced by mesocosm volumes (larger volumes were more stable) influence assembly processes of aquatic bacterial metacommunities along a press disturbance gradient. By combining path analysis and network approaches, we found mesocosm size categories had distinct relative influences of assembly process and environmental factors that determined spatiotemporal bacterial community composition, including dispersal and species sorting by conductivity. These processes depended on, but were not affected proportionately by, mesocosm size. Low fluctuation, large mesocosms primarily developed through the interplay of species sorting that became more important over time and transient priority effects as evidenced by more time-delayed associations. High fluctuation, small mesocosms had regular disruptions to species sorting and greater importance of ecological drift and dispersal limitation indicated by lower richness and higher taxa replacement. Together, these results emphasize that environmental fluctuations influence ecosystems over time and its impacts are modified by biotic properties intrinsic to ecosystem size.


INTRODUCTION
The community composition of both micro-and macro-organisms at a given point in space and time results from the interaction of multiple assembly processes, including ecological drift, species sorting (environmental filtering), dispersal, and speciation [1][2][3][4][5]. Most observational metacommunity studies, however, focus only on spatial snapshots without considering temporal dynamics of community assembly and association networks, or historical contingencies [2,6]. Hence, we still lack knowledge about the underlying mechanisms and regulating factors that temporal dynamics encompass.
When species sorting assembles communities, their composition tracks changes in environmental conditions that occur in time and space [2,7]. However, environmental tracking can be hindered or disrupted [8]. Such asynchrony can lead to historical contingencies by priority effects (e.g., [8,9], which can occur during early community formation or when communities re-assemble following perturbation. An important consequence of priority effects is that they impede or delay environmental tracking enacted by species sorting. Environmental changes may influence temporal community assembly processes and the strength of this can be regulated by ecosystem size (e.g., [6]. Studies have shown that microbial communities exposed to disturbances are initially, and often to a strong degree, stochastically assembled, but that the importance of species sorting increases later during community re-assembly as more species from the regional species pool arrive [10][11][12][13]. Rapidly fluctuating environmental conditions, however, may continuously disrupt environmental tracking by reducing opportunities for species sorting to select and shape local communities before the environmental conditions change again. This might promote coexistence of species with different niche optima [14][15][16] and, thus, reduce beta diversity [17], or could cause extinctions that bolster dispersal limitation and priority effects [18]. Nevertheless, many studies happen in controlled settings; thus, we lack knowledge on the temporal dynamics of these processes within larger, more complex habitats which track environmental changes [2,6]. Disturbance strength may uniquely affect microbial communities in ecosystems of different sizes as ecosystem size may influence assembly processes by increasing habitat heterogeneity, community abundance [6,19,20], and the pace at which communities track environmental changes. For instance, communities may experience different environmental variability including press disturbances (e.g., climate warming, eutrophication, or saltwater incursion), periodic and stochastic environmental fluctuations, where the latter may influence community assembly in response to the former over time and space.
Here, we implemented an experiment with freshwater bacterial metacommunities to test how different ecosystem size-induced environmental fluctuations influence the temporal dynamics of community assembly mechanisms. We collected a 64-day time series from mesocosms that gave bacterial communities time to experience natural environmental fluctuations. Specifically, we setup a natural experimental landscape with mesocosms containing identical lake water that differed in volume, which induced differences in environmental fluctuation intensity among the mesocosms. We created a press disturbance by applying a salinity gradient in which each mesocosm in a volume category had a different salinity as it has been shown that salinity affects bacterial communities in many ecosystems (e.g., [21][22][23][24]. The differences in the degree of environmental fluctuation among the size categories relied on larger water masses requiring more time to match changes in the surrounding air temperature combined with the increase in surface area-to-volume ratio that muted changes in the amplitude of salinity and other variable concentrations due to relatively less evaporation and precipitation. We hypothesized that the importance of species sorting would increase over time in local communities of larger mesocosms that experience relatively minor environmental fluctuations because their communities will have sufficient time for species selection in response to the initial salinity. Second, other environmental changes occurring in mesocosms would be slow in large mesocosms and this would allow time for taxa to be recruited from internal and external dispersal sources and to become active. We expected that species sorting related to salinity differences across communities, i.e., at the metacommunity scale, promotes recruitment of taxa best suited to the salinity. Last, we hypothesized that stochastic and/or dispersal-related assembly processes should be more important in small mesocosms where communities experience strong environmental fluctuations that continually disrupt environmental tracking. We combined quantitative path analysis methods that aim to estimate metacommunity processes with a network approach that identifies environmental tracking patterns through local and time-delayed co-occurrences to provide insights into temporal dynamics of microbial ecosystems [25].

Experimental set-up
Three different sizes (24.5, 70, or 200 L with surface areas 0.21 m 2 , 0.27 m 2 , and 0.39 m 2 , respectively) of hard-shell polyethylene mesocosms with an inverted conical frustum shape were arranged in a field beside Lake Erken (16 per size category) and filled with 0.1 mm filtered lake water from Lake Erken in Sweden (59°51'N 18°35'E) (water properties in Supplementary Information). Mesocosms were seeded with 1 L of sieved and mixed surface sediments collected from Lake Erken at~0.5 m water depth.
To induce species sorting with a press disturbance, a salinity gradient was created for each size category of mesocosms using nitrate-and phosphate-free sea salt (Red Sea Aquatics Ltd, Verneuil-sur-Avre, France). The gradient ranged from freshwater (0 ‰) to 6 ‰ with the salinity increasing by a 0.4% increment from one salinity level to the next (rationale for range in Supplementary Information). Mesocosm water surface area and volume were proportional such that air or rain dispersal was proportional across size classes. Mesocosm sediment was also a recruitment source [26][27][28]. Equal mesocosm bottom surface areas allowed for equal recruitment independent of fluctuation category.

Monitoring and sampling
Mesocosms were monitored on days 1, 2, and 4, and then every fourth day for 64 days from July to September 2016. Monitoring included depth profiles of conductivity (to measure salinity changes) and temperature, and depth-integrated pH, chlorophyll-a, and colored dissolved organic matter (CDOM) fluorescence (see Supplementary Information for details). Weather data from Svanberga, Sweden (0.87 km southwest of the site) included daily precipitation and hourly air temperature (Swedish Meteorological and Hydrological Institute). Every eighth day, water was collected for total organic carbon (TOC), total nitrogen (TN), and total phosphorus (TP) and analyzed using established methods [29].
Water samples for enumerating microorganism cells were collected simultaneously with bacterial community composition (below) and preserved with sterile formaldehyde to 2.5% [30]. Samples were stained with SYTO 13 Green Fluorescent Nucleic Acid Stain (ThermoFisher Scientific), counted (CyFlow Space flow cytometer, Partec, Münster, Germany) and analyzed using FlowingSoft software (Perttu Terho). Cell abundance was calculated as cells mL −1 while total community size was calculated as cell abundance (mL −1 ) multiplied by the entire mesocosm volume.

Bacterial community composition
Mesocosm water (0.5 L, depth-integrated) was collected on days 1, 2, 4, 8, and every 8 th day thereafter for 64 days to assess community composition through 16 S rRNA amplicon sequencing of reverse transcribed 16 S rRNA (cDNA) to specifically detect active members [31]. Ongoing immigration from external sources was characterized using sterile air and rain traps (n = 3, SA = 96 cm 2 ) interspersed beside the mesocosm units. Sheltered air traps contained 200 mL sterile filtered water whereas rain traps were uncovered. Every 8 days, trap samples were collected, pooled, and processed with mesocosm samples for 16 S rRNA gene amplicon sequencing. Both the water samples from the mesocosms (for RNA analysis), and the air and rain immigration samples (for DNA analysis) (details below), were collected and filtered onto 0.2 µm pore-size filters (47 mm Supor-200 filters, Pall Corporation, Hampshire, UK) until 5 min or 0.5 L volume was reached. Filters were flash-frozen in liquid nitrogen and stored at −80°C. DNA from initial lake water and sediment used in the experiment was sampled to learn initial communities and seed banks that could provide a species pool from which to identify sources of active bacterioplankton (RNA analysis) in the water column over the duration of the experiment.

Data processing
Sequencing resulted in 35.6 million paired reads from 609 demultiplexed samples including 12 extraction and PCR negatives. Primers were removed from sequences using cutadapt v 2.7 ref. [35]. The DADA2 pipeline [36] was used for sequence processing and taxonomy assignment of Amplicon Sequence Variants (ASVs) using the SILVA v. 138.1 reference database [37] (Supplementary Information, Table S1).
For beta diversity analyses, ASVs with counts less than 10 were removed and samples were subsampled to a minimum of 5028 reads, retaining 7983 unique ASVs. Samples not meeting the 5028 reads requirement were excluded (Table S2). Both alpha and beta diversity datasets represented 99% coverage (the probability that another individual collected from the original community has already been sampled) [38,39]. Raw sequences are available in the European Nucleotide Archive (study accession number PRJEB26595).
Fluctuation magnitudes among mesocosm sizes. For each environmental variable, fluctuations data were analyzed using the mean of absolute differences of mesocosms in a size category between one date and the previous sampling date. For variables with depth profiles (conductivity and temperature), the absolute difference at each depth was used to calculate the mean change per mesocosm. The environmental variables dataset is in the DiVA repository [42]. To determine if the magnitude of changes differed between mesocosm sizes over time, nonparametric tests for repeated measures with an ANOVA-type statistic (ATS) were used (R package and function nparLD, ref. [43]). Mesocosms were assessed using principal components analysis (PCA) of original and absolute changes of environmental variables (both log-transformed) and fit with environmental vectors (Fig. S2).
Community composition, diversity, and recruitment. Non-metric multidimensional scaling (NMDS) with Bray-Curtis dissimilarities was used to visualize bacterial community composition and environmental variables. Shannon index and Pielou's evenness were calculated and richness was estimated using the package "breakaway" [44]. Temporal beta diversity differences in each mesocosm were evaluated by comparing each community with the previous using Jaccard pairwise dissimilarity values. Dissimilarity was partitioned between taxa turnover (taxa replacement) and community nestedness (chronological subsets of taxa) using package "betapart" v1.5.2 ref. [45]. Variation from each partition captured by mesocosm size was compared using PERMANOVA tests [46] with function adonis and 999 permutations.
Recruitment was evaluated by pooling each mesocosm's active ASVs across days; ASVs present on day one were removed from the pool leaving those recruited during the experiment. Recruited ASVs were matched with their source seed bank(s) based on DNA from sediment, initial lake water, air, and rain. Sources for unmatched ASVs were considered unknown. For each mesocosm, the percent of recruited ASVs was calculated, split into each source, and examined across the salinity gradient using Pearson's correlations.
Path analysis. To detect drivers of metacommunity dynamics, a spatiotemporal path analysis was used [47]. This method calculates dissimilarity for all community pairs sampled over time and space and estimates, as individual paths on this beta-diversity measure, the influences of spatial distance (Δspace), temporal distance (Δtime), environmental distance (ΔEnvi), mean community size (<CommSize > , cell abundance multiplied by mesocosm volume), and absolute differences in community size (ΔCommSize) and taxa richness (ΔRich). A positive relationship between differences in community size and differences in richness could be explained by nestedness and would increase community dissimilarity [47]. Bray-Curtis dissimilarity was used for the community dissimilarity matrix (β Bray-Curtis ). A permutation-based approach adjusted with Benjamini-Hochberg procedure indicated path significance. Model fit was assessed with the standardized root mean square residual (SRMR). The analysis was run separately for each mesocosm size using the days required for network analysis (Supplementary Information), with the sem function in R package "lavaan" [48].
Network analysis. To uncover local and time-delayed microbial associations and the extrinsic effects of environmental variables on bacteria, extended local similarity analysis (eLSA) was applied [25]. Given our temporal data, this approach detects undirected associations (e.g., without time delays), and associations where the change of one factor (a taxon or environmental variable) chronologically leads or follows another factor. For a link between taxa and environmental variables, the association type (delayed or non-delayed) can indicate tracking that is time-lagged due to transient priority effects, or simultaneous through species sorting, respectively. Associations were determined for each mesocosm size category using eLSA wherein mesocosms within a size category were used as replicates (n = 16). Thus, for a given timepoint, each environmental variable, including salinity, from the same size category was normalized using the 'percentileZ' method and p mix which uses the determined theoretical P-value followed by permutation testing (n = 1000). Because of the within-and across-size variability of bacterial communities (e.g., significant differences in taxa richness), we selected and analyzed only the core bacterial groups for each mesocosm size to make it comparable. Hence, networks used the 50 most abundant ASVs from each size category. eLSA (v1.0.2) was run over eight sampling time points, allowing for local similarity (LS) correlations between samples taken eight days apart (d = 1). Only strong correlations with |LS | value ≥ 0.5 and Q ≤ 0.01 were considered and visualized in Cytoscape v3.8.2 [49]. Network characteristics were calculated using the Cytoscape plugin NetworkAnalyzer [50]. See Supplementary Information for details on sample selection, dominant ASV abundances, and statistics.

RESULTS
Environmental fluctuations in mesocosms Environmental variable fluctuations corresponded with mesocosm size and reflected rainfall and air temperature (Fig. S1, Table S3, Fig. S2). Size categories experienced significantly different conductivity and temperature fluctuations. After four days small and medium mesocosm conductivity fluctuated more than large mesocosms (Fig. S1, Table S3). Mean temperature fluctuation increased inversely with mesocosm size (Table S3). Mesocosm depth profiles showed stable conductivity, but temperature decreased with depth in medium and large mesocosms (Fig. S3).
Mesocosm sizes differed in nutrient concentrations and the absolute change of other environmental variables (chl-a, CDOM, pH, TN, TOC, TP and cell abundance, Table S3) and most pairwise comparisons showed that the degree of change differed significantly between sizes with the greatest changes in small mesocosms. Absolute changes between sampling dates and individual timepoints grouped according to size (Fig. S2). Measured nutrients and conductivity positively correlated with decreasing mesocosm sizes (environmental vector correlations, p < 0.05). For water temperature, sampling date was more influential than mesocosm size. Cell abundances per equal volume (cells mL −1 ) increased over time and were highest in small and medium mesocosms (ATS, p < 0.001, Fig. S4A). However, the total abundance of cells per an entire mesocosm was lower in small than medium and large mesocosms starting from day 24 (ATS, all p < 0.002, Fig. S4B).
Community composition, diversity, and recruitment Bacterial community composition shifted with time and initial conductivity in all mesocosm sizes ( Fig. S5-S7). Diversity indices (estimated richness, Pielou's evenness, Shannon index) did not differ on the first day (two-way ANOVA, all p > 0.05), but over time all three indices differed by mesocosm size (Fig. 1, repeated measures ANOVA, all overall p ≤ 0.001; pairwise Bonferroni adjusted). All sizes differed significantly in bacterial richness which was lowest in small and highest in large mesocosms (all p ≤ 0.001). The more evenly distributed ASV abundances in large mesocosms widened the separation in Shannon indices between large and small or medium mesocosms, indicating a greater presence of dominant and/or rare taxa in smaller mesocosms (Fig. 1). Mesocosm size explained some variability in beta diversity from turnover (F model = 13, R 2 = 0.04, p ≤ 0.001) with communities in small mesocosms experiencing higher turnover by taxa replacement than those in large mesocosms (Wilcoxon Test, W = 4276, Bonferroni p.adj. = 0.02, Fig. S8A). Statistically, mesocosm size did not explain variability in nestedness, although communities in large mesocosms trended towards greater nested species loss (Fig. S8B).
Recruited ASVs as a proportion of total unique ASVs, had weak negative Pearson's correlations with the salinity press disturbance in small and large mesocosms (r = −0.35 and −0.39, p < 0.005, respectively, Fig. S9). Less than 15 % of recruited ASVs in each mesocosm were attributed to a known source. In all mesocosm sizes, recruitment from water declined significantly with increasing salinity (small: r = −0.90, medium: r = −0.85, large: r = −0.89, all p < 0.001). Recruitment from sediment showed different patterns across salinity levels in small and large mesocosms: it decreased in small mesocosms and was unchanged in large mesocosms (r = −0.70, p = 0.002; r = 0.45, p = 0.08, respectively). Sediment was typically the largest recruitment source in the most saline mesocosms. Air and rain recruitment was related to salinity level only in the medium mesocosms where it was weakly positively correlated (air: r = 0.52, p = 0.04, rain: r = 0.58, p = 0.02). The fraction of 156 recruited ASVs from both rain and air sources was not correlated with salinity (all p ≥ 0.1).

Path analysis
Bacterial metacommunities in mesocosms of different sizes experienced disparate relative influences from species sorting by environmental variation, demographic stochasticity, and dispersal limitation (Fig. 2). The model fit for small mesocosms was roughly twice that of medium and large mesocosms (Fig. 2).
Species sorting (ΔEnvi) had the most influential direct effect on community dissimilarity (β bc ) (Fig. 2). This effect was strongest in small mesocosms and similar in medium and large mesocosms (sum of absolute standardized estimates 0.925, 0.773, and 0.766, respectively), but all sizes had significant environmental distance and community dissimilarity relationships (Tables S4-S6). Small mesocosms had five significant relationships between community dissimilarity and environmental variables (conductivity, temperature, chlorophyll-a, TOC, and TN); large mesocosms had three (conductivity, temperature, and chlorophyll-a), and medium mesocosms had only conductivity. Conductivity correlated most strongly with community dissimilarity of medium, followed by large and small mesocosms. Significant correlations between temporal (Δtime) or spatial (Δspace) distance and environmental (ΔEnvi) distance were positive and increased with mesocosm size. The indirect effect of time on community dissimilarity through species sorting was apparent with all measured variables except conductivity.
Demographic stochasticity was indicated by significant negative relationships between mean community size (<CommSize > ) and community dissimilarity in all mesocosm sizes (Fig. 2, Tables S4-S6). Small mesocosms had the strongest influence by demographic stochasticity. All mesocosm sizes had positive correlations between temporal distance and community dissimilarity indicating additional demographic stochasticity. Relationship strengths differed with size: temporal changes had the greatest influence in large, then small, then medium mesocosms.
Dispersal limitation shown as a positive correlation between geographic distance and community dissimilarity appeared only for small mesocosms (Fig. 2). Large mesocosms had a significant negative correlation between geographic distance and community dissimilarity but this was considered an artefact of the linear modelling framework [47] and negligible compared with the relationship between space and community dissimilarity via the environmental variation pathway.
The path analysis for medium and large mesocosms also suggested an effect of taxa nestedness whereby communities form as subsets of original communities over time or space (Tables S4-S6). First, differences in community richness (ΔRich) positively correlated with community dissimilarity. This relationship was strongest in large mesocosms. Second, a significant positive relationship occurred between differences in community size (ΔCommSize) and richness in medium and large mesocosms.

Association networks
Association networks of the 50 most abundant ASVs (members of Actinobacteriota, Bacteroidota, Cyanobacteria, Planctomycetota and Proteobacteria) differed among the three mesocosm sizes (Fig. 3, Table S7). The number of total edges and ASV nodes increased with mesocosm size, and the proportion of delayed (time-shifted) associations were higher in larger mesocosms (small: 25.9 %, medium: 43.5 %, large: 44.5 %) (Table S7). Small mesocosms had the most ASVs (n = 18) that were unassociated with environmental variables and bacterial abundance while medium and large mesocosms had only 5 and 8 ASVs, respectively (Table S7). Salinity had no associations with any 'core' ASVs at LS | correlations of ≥ 0.5.
Association networks were quantitatively compared by mesocosm size with commonly used topological characteristics. Negative associations, average number of neighbors, and network density (the proportion of possible edges that are associated with nodes) increased with mesocosm size (Table S7). In contrast, network heterogeneity (unevenness of the number of connections per node) and network centralization (the concentration of centrality among the nodes) decreased with ecosystem size. When considering only taxa associations, small mesocosms had the least centralized network with more taxa displaying similar numbers of links (Table S7).

DISCUSSION
Here we show how differences in environmental fluctuation strengths due to differences in ecosystem, i.e. mesocosm, size influenced the temporal dynamics of community assembly in response to a salinity press disturbance ( Fig. 4). First, species sorting was generally the most influential process for all mesocosms but there were differences in how species sorting operated among mesocosm sizes at the community (path analysis) and individual taxa levels (association network analysis). These evaluations indicated that under low environmental fluctuations, dominant ASV populations were effective trackers of environmental conditions. When ecosystem size-induced environmental fluctuations were strong (i.e., small mesocosms), environmental tracking was disrupted. Second, the salinity press disturbance altered community composition, especially under stable conditions (i.e., larger mesocosms), through the recruitment of taxa from seed banks (at high salinity sediment was the primary identified source). Third, stochasticity and dispersal-related assembly processes (e.g., dispersal limitation) were generally more important for communities of small ecosystems. Overall, our study aligns with previous findings that differences in ecosystem size have ecological consequences that influence community assembly processes [51][52][53], but here we identifed this effect to derive from the environmental fluctuations created by ecosystem size differences and corresponding differences in species sorting effects.
Salinity press disturbance enforces environmental tracking Differences in the magnitude of salinity press disturbances induced clear compositional shifts within and across mesocosms over time.
This was expected as we used salinity to induce species sorting because it is an environmental factor that causes clear taxonomic differences in aquatic bacterial communities [24,[54][55][56][57]. However, there were disparities in the role of salinity for structuring bacterial communities in each mesocosm size. The path analysis and network analysis results indicated that species sorting patterns differed across mesocosm sizes and were altered by time. The direct effects of significant environmental variables with unidirectional influences (i.e., salinity and temperature) on species sorting were most influential in medium and large, stabler mesocosms. However, when variables prone to feedbacks (i.e., nutrients, see below) were included into the total environmental effect on composition, species sorting was greatest in small, highly fluctuating mesocosms. In contrast, the indirect effect of time on composition via species sorting increased with mesocosm size and was driven primarily by changes in all environmental variables except salinity (which changed minimally within a mesocosm compared to the spatial salinity gradient). This temporal pattern generally agreed with the network results of the 50 most abundant bacteria which showed that they best tracked multiple environmental variables over time in medium and large mesocosms. Interestingly, almost all ASVs directly linked to environmental variables but populations of core groups of taxa did not oscillate strongly with temporal salinity changes (although weak associations, |LS | < 0.5, occurred).
The relationship between mesocosm size and environmental variables with the potential for feedbacks (i.e., nutrients) in the path analysis, suggests potential bottom-up effects. Although the path analysis portrays nutrients as effect variables, they are also modified by microorganisms. Likely due to the greater sediment: water ratio and the potential for salts to release sediment-bound nutrients through ion exchange, small mesocosms had greater densities of bacterial cell abundances and greater water nutrient and chlorophyll-a concentrations. Algal blooms were also observed, potentially increasing labile DOM resources via phytoplankton exudates. Indeed, resource availability and primary production may have played a major role in species sorting. In high-nutrient lakes, phytoplankton biomass and bacterial community composition were related [58]. These conditions could increase competition which hinders synchrony between abiotic variables and taxa [59].
There were clear taxonomic shifts in bacterial communities related to salinity preferences in space regardless of mesocosm size. For example, in agreement with other studies (e.g., [60]) Frankiales and Burkholderiales were relatively abundant only at low salinity levels whereas Synechococcales, which includes halotolerant species [61], became relatively more abundant at higher salinities (Fig. S7). However, the lack of a strong direct species sorting effect from temporal changes in salinity could result from several factors. First, salinity differed more across the salinity gradient (spatial changes) than within a mesocosm (temporal changes). Second, in the network analysis, we calculated associations only among the 50 most abundant taxa, thus, we likely overlooked conditionally rare taxa that can be temporarily abundant [62] as a consequence of the rapidly changing environment in small mesocosms. This is supported by the trend of higher taxa turnover and direct demographic stochasticity (discussed below) in small mesocosms. Last, bacterial communities can be an imprint of past environmental conditions [63] and the correlations detected between community dissimilarity and environmental variables might coincide with prior processes.
Taken together, our findings (conceptualized in Fig. 4) around the importance of species sorting and the strong temporal influences highlight the distinct differences in the mechanisms underlying species sorting in mesocosms of different sizes. These findings became apparent through combining the path analysis, which captures both spatial and temporal patterns at the whole community level, and the network analysis, which captures timeassociated patterns of the most abundant populations.
Ecosystem size indirectly regulates community assembly and associations among bacterioplankton While the different environmental conditions from the press disturbance and ongoing environmental changes throughout the experiment might explain why species sorting was the main driver of metacommunity assembly, our study suggests that other factors related to ecosystem size (e.g., spatial environmental heterogeneity within a mesocosm) could further regulate metacommunities.
The importance of species sorting can increase with environmental heterogeneity, i.e., the number of niches that are available for colonization [64,65]. In our study, large mesocosms contained greater within-mesocosm spatial environmental heterogeneity as evidenced by depth associated changes in temperature and light. This could explain why species sorting was more apparent (i.e., more associations between ASVs and abiotic variables) in large compared to small mesocosms in the network analysis. This increase may be attributed to (i) the greater availability of niches (and consistency of nutrients) found in larger mesocosms, or (ii) the synchronous establishment of bacteria which might have a better chance in a stable environment. In a study of protists experiencing light-dark fluctuations in aquatic microcosms and models, high fluctuations disrupted species synchrony between patches [66]. Within-mesocosm spatial heterogeneity could also explain the greater bacterial richness as ecosystem size increased.
Network topological features were partially influenced by mesocosm size: bacteria were more connected in medium or large than small mesocosms, suggesting that abundance dynamics were less similar in small mesocosms and indicating asynchrony among dominant bacteria. In the less densely populated, large mesocosms, competition may have been lower, which can lead to greater synchrony between species due to changes in abiotic conditions [59]. In our mesocosms, lower fluctuation and greater stability of the bigger mesocosms resulted in denser, more connected networks. Because size replicates in the network analysis spanned a salinity range that was wider than any range within a single mesocosm, the network analysis showed only weak tracking of the salinity changes over time by dominant bacteria. Specifically, salinity did not have strong synchrony with any core ASVs. This may result from methodological constraints where normalized salinity of withinsize mesocosms was used to assess associations. This indicates that the network approach could better detect environmental factors associated with species if they more equally (homogeneously) affected communities in space. Since studies often pool samples Fig. 4 Conceptual figure for the interpretation of statistical results and patterns based on path analysis, network analysis, and the partitioning of beta-diversity. In our study, the dominant deterministic force was the applied salinity press disturbance. Ecosystem size was manipulated by different volumes of mesocosms. Darker shading in bars indicates greater influence of the process.
over time or distant sampling sites, it can generate some biases and potentially mislead network inference.
Taken together, we suggest that these patterns indicate mesocosm size-specific mechanisms of species sorting: in small mesocosms, changes in community composition from species sorting primarily occurred through taxa replacement in response to variation in multiple environmental factors. In contrast, in larger mesocosms, environmental change was more gradual and cascaded into compositional differences through abundant bacteria tracking environmental changes over time by changing in relative population size, with lower replacement (Fig. 4).

Bacterioplankton recruitment
Initial community size differences due to mesocosm volumes might have affected subsequent community compositions through species sorting, but other factors related to the experimental set-up were unlikely to have substantial influence. The set-up ensured no extensive differences in the recruitment of novel species from external sources and estimated richness of active bacteria was equivalent on the first day. The dispersal sources (rain and air deposition, and seed banks in sediments and lake water) harbored high diversity and in previous studies were important recruitment sources for novel taxa following salinity disturbances [26][27][28] and other environmental changes [67]. Although large mesocosms contained more microorganisms and possibly a larger planktonic seed bank from which taxa could respond to the salinity disturbance, recruitment from water seed banks declined with salinity in all mesocosm sizes (Fig. 4). Even with reduced dispersal, future studies that extend beyond the 64 days sampled here may eventually see eco-evolutionary processes such as increased environmental tracking in small mesocosms due to bacterial diversification, which can be intensified by a history of environmental adversity [68]. The high percentage of ASVs with no identified source (85%) could indicate dispersal from other sources such as the snails we observed on most mesocosms or the effect of sequencing depth which can miss rare taxa. Additionally, sequencing depth could have incorrectly designated some taxa as recruited. We cannot discount that rare ASVs not detected on day one would be mistakenly categorized as recruited if detected later. However, this limitation should be consistent among treatments.

Roles of stochasticity and dispersal-related processes
Demographic stochasticity (leading to ecological drift) was an important driver of community assembly of all mesocosms via community size with the strongest direct effect in small mesocosms (Fig. 4). This result is bolstered by previous studies showing that ecological drift more often occurs in small communities [69,70] especially when the importance of species sorting is weak [71] or when the effective community size is small due to dispersal limitation [72]. This may be why we detected weak synchronous environmental tracking from the dominant populations across the small mesocosms. Drift can also alter the outcome of niche selection [73]. Nevertheless, the effect of time on community composition indicated that large mesocosm communities were most influenced by demographic stochasticity arising from temporal influences. In this case, large mesocosms may more strongly reflect (i) random changes in births and deaths from a community that grew in number over time, (ii) stochasticity based on priority effects from slower time-delayed tracking, or (iii) may reflect sampling timepoints that underrepresented the larger total community. Dispersal limitation as a driver of metacommunity dynamics (considering all mesocosms at one time point) was present only in small sized mesocosms and suggests that multiple communities emerged from similar initial conditions in the small mesocosms. However, the interpretation of the dispersal limitation is ambiguous (e.g., [74]). It could be true dispersal limitation whereby niche spaces that open (i.e., when species become inactive in response to the initial salinity changes which increase habitat specialists [75] and/or the strong environmental changes) remain empty [76]. However, it does not necessarily indicate true dispersal limitation between patches [74] or reduced immigration from a regional pool. Instead, it could be explained by low richness decreasing the likelihood that communities contain superb dispersers. When dispersal rates are low, local adaptations to environmental fluctuations can further strengthen priority effects by preemptive taxa [18]. In our study we found some evidence of transient priority effects indicated by a high number of time-lagged associations. However, they occurred only in larger mesocosms, whereas they might have been prevented by the overall higher taxa turnover in small mesocosms. Nevertheless, with our data and the applied approaches, it is not possible to clearly support or exclude dispersal limitation or priority effects and other factors that regulate them.

CONCLUSIONS
The novelty of our study is that we could show that the trajectories of (meta)community development are influenced by size-induced environmental fluctuations in concert with a salinity press disturbance. Overall, our results partially align with those from previous studies which show that after disturbances, stochastic community assembly initially is important, but the dominant influence shifts to deterministic processes in later successional stages (e.g., ref. [10]), especially when environmental conditions are stable. Dispersal limitation and ecological drift (demographic stochasticity) were drivers of metacommunity dynamics after community establishment with strong environmental fluctuations. Moreover, our results indicate that mesocosms with reduced environmental fluctuations may facilitate considerable time-delayed species sorting, potentially due to transient priority effects. Collectively, our study highlights that environmental fluctuations, resulting from the temporal environmental change dynamics, are important to consider in future community assembly studies. Future studies should also include top-down effects that are potentially altered by ecosystem sizeinduced environmental fluctuations in addition to bottom-up effects that we focused on here.
The importance to account for the interactive effects of environmental fluctuations and ecosystem size on community assembly can also have implications for different ecosystems that are affected by both the magnitude of environmental fluctuations and additional size-dependent properties of their ecosystem, such as different types of small water bodies, animal guts, or bioreactors. Further, the framework presented here may inform predictions of how drought or irrigation that affect aquatic ecosystems can potentially reduce species sorting synchrony and enhance maladapted taxa.

DATA AVAILABILITY
The data supporting the results are archived in the public repository European Nucleotide Archive with accession number PRJEB26595 and environmental data are made available in the Swedish institutional repository, DiVA, (diva-portal.org) with the following accession number: diva2:1210995.