Spatiotemporal assembly and functional composition of planktonic microeukaryotic communities along productivity gradients in a subtropical lake

Microeukaryotes play crucial roles in the microbial loop of freshwater ecosystems, functioning both as primary producers and bacterivorous consumers. However, understanding the assembly of microeukaryotic communities and their functional composition in freshwater lake ecosystems across diverse environmental gradients remains limited. Here, we utilized amplicon sequencing of 18S rRNA gene and multivariate statistical analyses to examine the spatiotemporal and biogeographical patterns of microeukaryotes in water columns (at depths of 0.5, 5, and 10 m) within a subtropical lake in eastern China, covering a 40 km distance during spring and autumn of 2022. Our results revealed that complex and diverse microeukaryotic communities were dominated by Chlorophyta (mainly Chlorophyceae), Fungi, Alveolata, Stramenopiles, and Cryptophyta lineages. Species richness was higher in autumn than in spring, forming significant hump-shaped relationships with chlorophyll a concentration (Chl-a, an indicator of phytoplankton biomass). Microeukaryotic communities exhibited significant seasonality and distance-decay patterns. By contrast, the effect of vertical depth was negligible. Stochastic processes mainly influenced the assembly of microeukaryotic communities, explaining 63, 67, and 55% of community variation for spring, autumn, and both seasons combined, respectively. Trait-based functional analysis revealed the prevalence of heterotrophic and phototrophic microeukaryotic plankton with a trade-off along N:P ratio, Chl-a, and dissolved oxygen (DO) gradients. Similarly, the mixotrophic proportions were significantly and positively correlated with Chl-a and DO concentrations. Overall, our findings may provide useful insights into the assembly patterns of microeukaryotes in lake ecosystem and how their functions respond to environmental changes.


Introduction
Microeukaryotic plankton, encompassing a broad range of unicellular eukaryotes like algae, protozoa, and fungi, are vital components within microbial food webs (Azam and Malfatti, 2007).More importantly, the microplankton community contributes significantly to primary productivity and nutrient recycling in freshwater environments (Liu et al., 2013;Gad et al., 2020).The functions and community dynamics of planktonic microeukaryotes exert a profound influence on water quality (Yu et al., 2014;Gad et al., 2022).Adverse events like algal blooms and toxin production, induced by community variation, adversely affect sources of drinking water (Zhang et al., 2018;Treuer et al., 2021).Therefore, planktonic microeukaryote communities are commonly employed as biomonitors for assessing aquatic ecosystem health and water quality (Foissner and Berger, 1996).Over the past decades, much effort has been made to unveil how the freshwater microeukaryotic communities are distributed and assembled (Chen et al., 2019;Isabwe et al., 2019;Abdullah Al et al., 2022).The biogeography of microeukaryotic planktons can be shaped by abiotic factors (e.g., nutrients, temperature, light, pH, depth, spatial factors) and biotic factors (including top-down predation, competition, mutualism, and trade-offs) (Gong et al., 2015;Wang et al., 2015;Chen et al., 2019;Wang et al., 2020).Nevertheless, discerning the relative significance of these factors in determining community structure and understanding the variability of those driving factors across spatial, vertical, and seasonal scales, remains fundamental issues in freshwater ecology.Furthermore, understanding how freshwater microeukaryotic community dynamics evolve is critical for ecosystem functionality prediction and management in the context of global change.
Pigmented microeukaryotes are key constituents of the phytoplankton realm in aquatic ecosystems (Xu et al., 2018).The concentration of chlorophyll a (Chl-a), a reliable surrogate for photosynthetic potential and primary productivity, often demonstrates a hump-shaped relationship with the diversity of phytoplankton in aquatic environments (Vallina et al., 2014).The productivity-diversity relationships depend on scale of investigations (Chase and Leibold, 2002), the type of water body (Borics et al., 2014), and are susceptible to the influence of environmental change (Righetti et al., 2019).Studies have documented that phytoplankton release a substantial proportion of photosynthate, ranging from 0 to 80% of that fixed via photosynthesis, as dissolved organic matter (DOM) into the surrounding waters (Hulatt and Thomas, 2010).The DOM is directly incorporated into the microbial loop via heterotrophic bacterial consumption, subsequently enhancing the abundance of bacterivorous (i.e., heterotrophic/mixotrophic) microeukaryotes.In a recent investigation, an uncommon U-shaped pattern was discovered between microeukaryotes and Chl-a in eutrophic coastal oceans (Wang et al., 2020).Nevertheless, within freshwater ecosystems, there is a lack of comprehensive exploration into the variations in microeukaryotic community structure linked to productivity and the primary environmental factors responsible for driving these changes.
Identifying and quantifying the ecological processes that regulate the assembly of microbial communities in aquatic ecosystems stands as a central focus in microbial ecology (Logares et al., 2013;Nemergut et al., 2013).Addressing this critical issue is essential for advancing our comprehension of diversity, functioning, and successional dynamics of microorganisms (Liu et al., 2019).Classical niche-based theory proposes that deterministic processes, involving both biotic interactions and abiotic factors, jointly influence community structure (Chesson, 2000;Fargione et al., 2003).However, in neutral theory of biodiversity, species are viewed as functionally equivalent, with their abundances primarily determined by stochastic processes such as birth, death, speciation, extinction, and dispersal (Hubbell, 2001;Chave, 2004).Deterministic processes mainly governed the assembly of microeukaryotic communities in estuarine, coastal, and marine ecosystems (Wu et al., 2018;Zhu et al., 2021;Wu et al., 2023).However, in freshwater environments such as rivers and lakes, as well as surface ocean, microeukaryotic structure exhibited a pronounced influence of stochastic processes (Chen et al., 2019;Logares et al., 2020;Ren et al., 2022;Han et al., 2023).Frequently, microeukaryotic communities seem to be regulated by a combination of both deterministic and stochastic processes (Gad et al., 2020;Hou et al., 2020;Mo et al., 2021).Nonetheless, a comprehensive understanding of microeukaryotic geographical patterns and assembly mechanisms in freshwater ecosystems across spatial, temporal, vertical water column, and environmental gradients is still lacking.
The Qianxia Lake (E119.9N28.1), located in Qingtian County, Zhejiang Province, China, is a fjord-type lake with a long and narrow watershed covering approximately 3,330 km 2 .It spans 105 km in length and averages 31.7 km in width (Yao et al., 2022).The lake is a crucial source of drinking water for numerous towns and sustains agriculture, power generation, and tourism, with minimal industrial pollution.Thus, the Qianxia Lake watershed offers an ideal setting to explore the relative importance of environmental and spatial factors on community composition.This knowledge is crucial for integrating planktonic microeukaryotic components into ecosystem and biogeographical models to improve predictive accuracy.Furthermore, given the dispersal capacity of microplankton communities in lotic habitats, it renders them well-suited for investigating the impact of dispersal on community variation and the influence of seasonal changes on dispersal.
In this investigation, we conducted a comparative analysis of microeukaryotic plankton communities in the Qianxia lake, employing amplicon sequencing of 18S rRNA genes.Our data included 60 water samples collected across two seasons and various water depths along the Qianxia Lake (Supplementary Figure S1), with the following objectives: (1) Unraveling the spatiotemporal distribution of microeukaryotic communities and deciphering the underlying mechanism governing their assembly; (2) Investigating how the microeukaryotic diversity and community structure respond to productivity gradients; and (3) Assessing the impact of environmental factors on variation in functional composition of microeukaryotic communities.

Sampling and determination of environmental parameters
Sampling of water was carried out at 8 sites in spring and 12 sites in autumn across the Qianxia Lake (119°65′-120°04′ E, 27°99′-28°13′ N) in April and September 2022 (Supplementary Figure S1).A total of 60 water samples were collected at depths of 0.5, 5, and 10 m in the epilimnion zone of water column.Due to the lake's depth ranging from 10.5 to 79.4 m, deep samples were collected away from the lakebed to prevent disturbance to sediment (Supplementary Table S1).The water collection utilized 5 L Niskin bottles, pre-filtered through a 200 μm mesh to eliminate larger plankton and debris.Subsequently, water (~500 mL) with microeukaryotes (size <200 μm) was vacuumfiltered onto a 0.22 μm filter (Supor 200, PALL, Michigan, United States).All membranes were placed into 2 mL cryotubes and preserved in liquid nitrogen for DNA extraction.
The in situ lake water's physicochemical parameters (temperature, pH, and dissolved oxygen) were determined through a water quality analyzer (HQ40d, Hach company, United States).The Chl-a concentration was measured employing an electronic sensor (TriBox mini, Germany).Determination of nitrate (NO 3 -N), nitrite (NO 2 -N), ammonium (NH 4 + -N), total nitrogen (TN), and soluble reactive phosphate (PO 4 3− ) in all sub-samples via an auto-analyzer (Seal, Germany).Total phosphorus (TP), COD Mn , and total organic carbon (TOC) were determined through spectrophotometry following standard methods.A total of 12 environmental parameters were determined in this study (Supplementary Table S1).

Bioinformatics
The raw sequence data were initially demultiplexed, and primers were trimmed with Cutadapt v1.18 (Martin, 2011), followed by processing through the DADA2 pipeline (Callahan et al., 2016) in R v4.0.0 (R Core Team, 2018).The parameter settings followed our prior study (Zou S. et al., 2021).Briefly, sequencing reads were quality filtering with the "filterAndTrim" function.The dereplication process generated unique sequences, implemented with "derepFastq" function; Error models were trained with the "learnErrors" function.ASVs were inferred using the core function "dada." Paired reads were merged using the "mergePairs" function, and bimeras were eliminated with the "removeBimeraDenovo" function.Taxonomic classification was performed using the DECIPHER package's "IdTaxa" function against the PR 2 v2.0 database (Wright, 2016).ASVs that were unassigned or appeared as singletons, as well as reads from Metazoa, Ulvophyceae, Streptophyta, and Rhodophyta, were excluded from subsequent analysis.
Prior to computing alpha diversity metrics, all sample were rarefied to the lowest sequencing depth of 18,300 reads (Supplementary Table S3), which minimize sequencing depth bias and enable comparable diversity assessments.To analyze beta diversity analysis, we normalized the ASV table using edgeR package (Robinson et al., 2010), following the recommendation of McMurdie and Holmes (2014).Bray-Curtis dissimilarities distances were computed and then visualized in NMDS plot using the "metaMDS" function of the vegan package (Oksanen et al., 2013).

The relative contribution of environmental and spatial factors in structuring microeukaryotic communities
We assessed the spatial distribution of microeukaryotic communities by assessing the community turnover rate (z value) over spatial distances, as proposed by Ranjard et al. (2013): log 10 (χd) = (−2z) * log 10 (d) + b.The pairwise Sørensen similarity (χd) was computed using the "dsvdis" function from the labdsv package (Roberts, 2019).The distance (d) between sampling sites (in meters) was determined using the function "distm" of geosphere package (Hijmans et al., 2017).The "b" represents the intercept of the linear model.
We further employed variation partitioning analysis (VPA) to assess the contribution of environmental and spatial factors in microeukaryotic assembly (Borcard et al., 1992).All environmental factors (except pH) underwent square root transformation for normality and homoscedasticity.For spatial factors, we initially converted the longitude and latitude to the Cartesian coordinates using the function "geoXY" in the SoDA package (Chambers, 2008).Subsequently, an Euclidean distance matrix was generated from these Cartesian coordinates using the "dist" function.The "PCNM" function (permutations = 1,000, vegan package) was then applied to this matrix (Borcard and Legendre, 2002).
Redundancy analysis (RDA) was chosen to explore the relationship between microeukaryotic communities and environmental/spatial factors based on the longest gradient length of detrended correspondence analysis (DCA).Preceding RDA, the following steps were executed: (1) Hellinger transformation of microeukaryotic data; (2) Variance inflation factor (VIF) was calculated using the "vif.cca"function (vegan package) and variables with VIF > 10 were discarded to address the multicollinearity effects; and (3) Identification of significant (p < 0.05) explanatory factors through forward selection using the "ordiR2step" function (Blanchet et al., 2008).The relative effects of environmental and spatial factors on community variation were determined with VPA using the "varpart" function of package vegan in R, and the significance of each fraction was determined using the permutation test with the "anova.cca"function.

Functional analysis based on taxonomic specifically traits approaches
To elucidate the functional composition of microeukaryotes and their responses to the environmental gradients, we utilized a trait-based approach, as previously described (Genitsaris et al., 2015;Adl et al., 2019;Ramond et al., 2019;Wang et al., 2020;Pan et al., 2022).In brief, the taxonomic assignments of microeukaryotic members, as annotated by the PR 2 database, were converted into ecological traits, while retaining their read proportions within a given community.In this study, the  S4).The ASV richness and abundance for these three functional groups were determined by summing the taxonomic groups within each corresponding functional category.

Estimation of importance of environmental variables using random forest model
We utilized random forest analysis to evaluate the impact of environmental factors on the functional composition of microeukaryotes, employing the rfPermute package with 999 permutations (Archer, 2019).The analysis of percentage increases in mean squared error (MSE%) of predictors determined the importance of the variables.Variables with higher MSE% values were implied as more influential.Model significance and cross-validated R 2 were assessed using the A3 package in R (Fortmann-Roe, 2015).

Fitting the neutral community model for microeukaryotes
We evaluated the influence of stochastic processes on the microeukaryotic community assembly using the neutral community model (NCM), as described by Sloan et al. (2006).The NCM employs R 2 for the fit to the neutral model and Nm values as the product of metacommunity size (N) and immigration rate (m).Additionally, we compared the richness and abundance of microeukaryotes in neutral and non-neutral (above and below) partitions to evaluate the deviations from the NCM predictions.All computations and visualizations were performed using R.

Statistical analyses
Two-tailed student's t-tests and one-way ANOVA with least significant difference (LSD) post-hoc were utilized to assess the variations in physicochemical variables, alpha diversity metrics, relative proportions of major taxonomic taxa between seasons, among depths, and across Chl-a levels.Permutational multivariate analysis of variance (PERMANOVA, 999 permutations) was employed to examine differences in community structure among groups based on Bray-Curtis distances (Anderson, 2001).All geographic measurements were calculated and visualized with ArcGIS v.10.8 (Redlands, California, United States).The relationships between environmental or spatial factors and taxonomic/functional diversity were assessed through Pearson or Spearman's rank correlations analysis using SPSS v.13.0 (IBM, United States).

Differences in alpha diversity and community composition of planktonic microeukaryotes across seasons and depths
A total of 6,383 ASVs were generated from 6,520,323 high-quality sequences for 60 samples, among which 2,758 and 4,059 ASVs were detected in spring and autumn, respectively (Supplementary Table S3).Of the rarefied dataset, the alpha-diversity metrics of microeukaryotes exhibited considerable variability among samples, with ASV richness ranging from 107.3 to 850.9 and Shannon indices between 3.4 to 7.9 (Supplementary Table S3).In contrasting the two seasons, autumn had higher species richness (338.1)compared to spring (269.4), while Shannon indices lower in autumn than in spring.Among the three depths, the alpha-diversity estimators usually increased with the depths (Figures 1A,C).The species richness and Shannon both formed significant hump-shaped relationships with concentrations of Chl-a (R 2 ≥ 0.115, p ≤ 0.033), and the species richness appears to peak at the concentrations of Chl-a of 4.0 μgL −1 (Figures 1B,D).

Changes in taxonomic composition between seasons, among depths, and along productivity gradients
The NMDS plot and PERMANOVA analysis indicated a distinct separation in microeukaryotic community compositions between spring and autumn (R 2 = 0.302, p = 0.001).However, samples from different depths did not exhibit clear differentiation (R 2 = 0.025, p = 0.79; Table 1; Figure 3A).The log-transformed Sørensen similarity exhibited significant distance-decay patterns in microeukaryotic community for both seasons (adj.R 2 ≥ 0.203; p < 0.01, Figure 3B).The turnover rate (z value) was slightly higher in autumn than in spring (0.03 vs. 0.02, Figure 3C).Similarly, microeukaryotic community similarities between any two sampling sites significantly decreased with pairwise differences in Chl-a concentrations (adj.R 2 ≥ 0.299; p < 0.001 Figure 3D), highlighting the crucial role of Chl-a in structuring microeukaryotic assembly.Additionally, similarity of planktonic microeukaryotes, spanning taxonomic ranks from species to phylum, also demonstrated pronounced distance-decay patterns (p < 0.001; Supplementary Figure S3).
To explore the impact of the productivity gradients on microeukaryotic community structure, we statistically compared the relative proportions of major taxa at three Chl-a levels (0-2 μg L −1 for low level; 2-5 μg L −1 for intermediate level; and 5-11 μg L −1 for high level; Figure 5; Supplementary Table S5).Proportions of Cryptophyta, Dinoflagellata, Basidiomycota, and Apicomplexa increased with levels of Chl-a.Under intermediate levels, six taxonomic groups (Haptista, Ochrophyta, Ciliophora, Ascomycota, Choanoflagellida, and Tubulinea) had the highest relative abundance.However, Chytridiomycota and Bacillariophyta were more abundant in low concentrations of Chl-a (Figure 5).

Variations in functional structure of microeukaryotic communities
The ASV richness of phototrophs varied greatly across all samples (51-322) and was comparable to that of heterotrophs (41-386).Phototrophs and heterotrophs were both nearly three times more abundant than mixotrophs (Figure 6A).Seasonally, phototrophs and heterotrophs showed similar patterns, with significantly higher abundances in autumn than in spring (p ≤ 0.04).Conversely, mixotrophs richness was slightly higher in spring, although this difference did not reach statistical significance (p = 0.648; Figure 6B).Regarding depth-related comparisons, a notable difference in richness  Differences in diversity and composition of planktonic microeukaryotes community.(A) Pie plots depict ASV numbers for major taxa during spring and autumn, respectively.(B) Bar plot illustrates relative proportions of major taxa across the samples.was evident solely for heterotrophs between bottom layer and middle layer (55.5 vs. 49.4;p < 0.05; Figure 6B).The correlations between the richness of three functional groups and Chl-a concentrations were weak and lacked statistically support, except for mixotrophs (R 2 = 0.218, p < 0.001; Figure 6C).Seasonal variations significantly influenced the relative abundances of the three functional groups of microeukaryotes (Figure 7A).Heterotrophs and phototrophs dominated the microeukaryotic communities among all samples, accounting for an average percentage of 46.9 and 42.8%, respectively, while mixotrophs consistently had a low contribution (10.2%).Higher abundances of mixotrophs and heterotrophs were detected in spring, while phototrophs became more abundant in autumn (p = 0.001).However, vertical comparisons did not reveal any significant differences in reads proportions of functional groups (p > 0.436; Figure 7B).
Random forest analysis was utilized to assess the significance of environmental factors in shaping functional composition.The results revealed the highest contribution to heterotrophs (36.6%), followed by mixotrophs (30.5%), and phototrophs (23.9%).These variations were predominantly driven by a combination of environmental factors, including Chl-a, N:P ratio, DO, TP, TOC, temperature, and NO 2 -N (p < 0.05, Figure 8A).Relationships between the relative abundances of functional groups and key environmental factors exhibited distinct trends: the proportion of phototrophs showed a significant and positive correlation with N:P ratio (R 2 = 0.097, p = 0.016), whereas heterotrophs exhibited a significant and negative correlation with N:P ratio (R 2 = 0.157, p = 0.002).Nevertheless, the correlation between mixotrophs and N:P ratio was weak and non-significant (R 2 = 0.031, p = 0.184; Figure 8B).Across the gradients of Chl-a, phototrophs, heterotrophs, and mixotrophs exhibited hump-shaped, U-shaped, and positive linear patterns, respectively.Similarly, these relationships were in parallel with dissolved oxygen gradients but with lower correlation coefficients, and only significant for phototrophs and mixotrophs (Figures 8C,D).

Relative contribution of environmental drivers to variations in microeukaryotic communities and underlying assembly mechanisms
RDA demonstrated variations in the percentage of explained variance for microeukaryotic communities by the first two axes across FIGURE 3 nMDS using Bray-Curtis dissimilarity illustrates the seasonal variations in planktonic microeukaryotic communities (A).Scatter plots depict distancedecay patterns between Sørensen similarity of communities and pair-wise geographic distance in spring, autumn, and both seasons (B,C).Sørensen similarities of communities also show a negative correlation with the gradients of productivity (D).Zou et al. 10.3389/fmicb.2024.1351772Frontiers in Microbiology 08 frontiersin.orgseasons, amounting to 39.5, 41.96, and 41.18% for spring, autumn, and both seasons, respectively (Figures 9A,B; Supplementary Figure S4).Environmental (Chl-a and NO 2 -N) and spatial variables (PCNM1, PCNM2, and PCNM5) exerted notable impacts on microeukaryotic community composition in both seasons by forward selection.In autumn, pH, PCNM3, PCNM4, and PCNM7 also showed significant correlations with communities.Variation partitioning analysis (VPA, Figures 9C,D) indicated that purely spatial variation played a larger role in community composition during spring (14.9%), autumn (21.6%), and both seasons (34.3%) compared to purely environmental The relative proportions of microeukaryotic taxa with significant differences between spring and autumn (A-G) or among vertical depths of water column (H-J).Significance was assessed through Student's t-test or one-way ANOVA test (p < 0.05).The error bars denote standard errors.
factors (0.9% in spring, 2.5% in autumn, 9.3% in both seasons).Additionally, more than half of the proportion of community variation for each and both seasons was unexplained (51.3-64.2%,Figures 9C,D; Supplementary Figure S4).
To assess the influence of stochastic processes on the assembly of microeukaryotic communities during spring and autumn, we employed the neutral community model (NCM).Results demonstrated that the NCM explained a substantial fraction (62.8, 66.7, and 55.1%) of variation in community composition for spring, autumn, and both seasons combined, respectively (Figure 10).This underscores the critical role of stochastic processes in influencing the assembly of microeukaryotic communities.The N m value for microeukaryotic taxa was lower in spring (N m = 102) than in autumn (N m = 337), indicating that higher species dispersal in autumn (m: 1.84% vs. 0.56%), due to sequence depth that was rarefied for both seasons (N = 18,300).Furthermore, the proportions of microeukaryotic Differences in community richness among samples for phototrophic, mixotrophic, and heterotrophic groups (A).Comparisons between two seasons and among vertical depths (B).The mixotrophs formed a significant humped relationship with chlorophyll-a concentrations, contrasting with the heterotrophs and phototrophs (C).T-test and ANOVA were employed to assess differences in richness between seasons and among depths.Photo, phototrophs;Mixo, mixotrophs;Hetero, heterotrophs. 10.3389/fmicb.2024.1351772Frontiers in Microbiology 10 frontiersin.org ASVs revealed that the neutral fraction contributed significantly more to species richness (88.7-94.6%)than the above and below fractions in both seasons, but not to ASVs' abundance (17-54.9%; Figure 10D).The findings, along with VPA results, suggest that stochastic processes primarily influenced community structure.

Discussion
Strong seasonality of planktonic microeukaryotic communities with negligible depth influence The nMDS plot distinctly revealed the separation of samples between spring and autumn, indicating significant seasonality in planktonic microeukaryotic composition (p < 0.05) (Table 1; Figure 3).The seasonal differentiation can be attributed to variations in physicochemical factors, biotic interplays, and hydrologic environment (Sommer et al., 2012).Firstly, we observed significant differences in specific environmental factors (temperature, DO, pH) and nutrients (NO 2 -N, PO 4 -P, TP, and TOC) between the two seasons (Supplementary Table S2).Water temperature, a key seasonal factor, may directly regulate microeukaryotic community composition (Liu et al., 2013).Prior investigations indicated that elevated water temperatures could enhance the proliferation of microeukaryotic species (Liu et al., 2020), partially explaining the higher species richness in the autumn compared to spring (Figure 1A).Nitrite (NO 2 -N) emerged as the sole nutrient significantly co-varied with microeukaryotic community structure during both seasons, as revealed by RDA ordination (Figures 9A,B).Nitrite is a vital form of inorganic nitrogen readily assimilated by phytoplankton.Certain phytoplankton populations have been reported to consider nitrite as a more important nitrogen source than nitrate (Yves, 1998).Hence, nitrite concentration could mediate the phytoplankton community variation, thereby indirectly influencing the entire microeukaryotic community composition within the lake system.Secondly, the seasonal dynamics of microbial plankton community represent a repeatedly recurring process primarily influenced by environmental fluctuations and internal interactions, such as competition and predation (Sommer et al., 2012).Thirdly, the strong seasonality could also be explained by the high turnover of microeukaryotic communities (Zhang et al., 2017), with stochastic processes (ecological drift and dispersal) playing a significant role in structuring their communities (Logares et al., 2018).The spring-autumn variations in microeukaryotic community composition were mainly characterized by a higher abundance of Dinoflagellata, Basidiomycota, Mesomycetozoa, and Apicomplexa in spring, while Chlorophyta, Haptista, and Choanoflagellida were more abundant in autumn (Figure 4; Supplementary Table S5).It is worth noting that the protists, such as ciliates, dinoflagellates, and fungi, possess much higher SSU rDNA copies than other protists, which could potentially result in an overestimation of their prevalence in environmental surveys (Godhe et al., 2008;Gong et al., 2013).Variations in single-celled rDNA copy number have been tightly linked with cell size, thus it is more appropriate to interpret the rDNA reads in relation to relative proportion of eukaryotic biomass (Fu and Gong, 2017;Zou S. et al., 2021).Therefore, planktonic dinoflagellates and fungi (Basidiomycota and Mesomycetozoa) might bloom in spring, resulting in higher biomass during this season.
Limited variability was noted in planktonic microeukaryotic communities across water column depths (Figure 3A; Table 1), despite there being significant differences in certain environmental variables (e.g., DO and temperature; Supplementary Table S2).This observation aligns with a previous study on bacterioplankton communities in a reservoir with a depth range of 0.5-20 m (Zlatković et al., 2022) but contrasts with research covering a much wider depth range.For example, Gong et al. (2015) documented that structure of microeukaryotic community was significantly shaped by depth in sediments of coastal shallow water (15-75 m).In addition, David et al. (2021) investigated the protistan community in Lake Baikal, Russia, and demonstrated a strong significant effect of depth (5-1,400 m) on protist community stratification.The limited variation in planktonic microbial communities within the epilimnion zone (0-20 m) of water column is likely due to wind-induced water mixing, which reduce stratification and constrain niche differentiation opportunities for microorganisms.Therefore, the homogeneity of microeukaryotic plankton can be understood as breakdown of stratification in the upper layer of water column, may have the following ecological consequences: (1) Upbringing nutrient-rich deeper waters to the surface, increasing the nutrient availability for diverse microbial groups and homogenizing environmental conditions (e.g., temperature, pH);  Species richness-productivity relationship in lake microeukaryotic communities follows a significant hump-shaped pattern The species richness and productivity relationships (SRPR) at regional and local scales are long-standing and fundamental aspects of ecological research.SRPR can exhibit various patterns, such as linear, unimodal, U-shaped, or non-significant patterns (Waide et al., 1999;Mittelbach et al., 2001;Whittaker and Heegaard, 2003;Wang et al., 2020).Among these patterns, hump-shaped relationships occur most frequently, particularly in aquatic systems (Mittelbach et al., 2001).Our current study demonstrated a hump-shaped curve between microeukaryotic diversity and Chl-a (Figures 1B,D), and this pattern was mainly attributed to Dinoflagellata, Basidiomycota, Oomycota, and Perkinsozoa, most of which were heterotrophs, with the exception of Dinoflagellata (Supplementary Figure S2; Supplementary Table S4).Our observations of microeukaryotic species richness maximizing at intermediate productivity levels (~4.0 μgL −1 ) align with earlier research in marine (Spatharis et al., 2008) and lake ecosystems (Dodson et al., 2000;Grover and Chrzanowski, 2004).Several hypotheses have been proposed based on existing ecological theories to explain the hump-shaped relationship.At the low productivity zone (left end of the curve), species richness seems primarily influenced by nutrient limitation (Spatharis et al., 2007).As productivity increases, the abundance of rare species may surge, leading to increased species richness (Abrams, 1995).Nonetheless, within the high productivity regime (rightmost of the curve), a reduction in species richness may be due to increased competition or severe environmental stressors (Dodson et al., 2000;Spatharis et al., 2008).These findings underscore the complex ecological dynamics that govern microeukaryotic biodiversity patterns in aquatic ecosystems.

Linking functional shifts of microeukaryotic planktons to environmental drivers
The diversity and taxonomic composition of communities exhibited significant correlations and formed hump-shaped relationships with Chl-a and DO (Figures 1B,D, 3D; Supplementary Figure S2).However, no similar patterns were observed for phototrophs and heterotrophs, except for mixotrophs when using trait-based functional analysis (Figure 6C).The Random forest analysis indicated that the alterations in functional composition of freshwater microeukaryotic communities result from various environmental factors (i.e., N:P ratio, Chl-a, DO, TP, and TOC; Figure 8), which are broadly consistent with findings from several recent investigations (Amorim and Moura, 2021;Lu et al., 2023;Xue et al., 2023).The determined nutrient levels showed relatively high contents of PO 4 -P (on average 0.14 mg/L) but low levels of DIN (dominated by NO 3 -N, 0.56 mg/L) in the regime of Qianxia Lake.Furthermore, the mean N:P ratio of 4.6:1 was significantly below the Redfield ratio of 16:1.Collectively, these results indicate an N-deficit condition within the studied region.The N limitation may be one of the reasons for both phototrophic and mixotrophic abundance positively correlated with the N:P ratio (Figure 8B).Under N-deficit condition, the elevated prevalence of heterotrophs could be explained by heterotrophic microeukaryotes obtaining nutrients by grazing bacteria (Sherr and Sherr, 2002).The stronger relationship of relative proportion with Chl-a rather than DO in this study suggests that productivity probably plays a more crucial role than DO in indicating the abundance of mixotrophs.This finding, however, contrasts with our previous investigation of pico−/nanoeukaryotes in coastal oceans (Wang et al., 2020).

Significant distance-decay relationship suggests pronounced dispersal limitation effects
In the current investigation, a notable distance-decay pattern of the microeukaryotic community was observed in Qianxia Lake of Southeast China (Figures 3B,C).The calculated turnover rate (z value) for microeukaryotes in our research was 0.014, falling within the range (0.0019 ~ 0.26) documented for microbial communities by Woodcock et al. (2006).This pronounced turnover in biodiversity is typically ascribed to a dynamic balance between extinction and immigration, coupled with an increase in habitat diversity over greater spatial distances (Kallimanis et al., 2008).Previous microbial biogeographic investigations have disclosed a decline within certain populations over distance, highlighting the substantial influence of dispersal limitation and resulting in conspicuous regional uniqueness (Martiny et al., 2006;Isabwe et al., 2018;Bai et al., 2022).This challenges the microbial ubiquitous dispersal hypothesis (Finlay and Clarke, 1999).Given the potential impact of both deterministic and stochastic factors on dispersion, dispersal limitation could theoretically be influenced by stochastic, deterministic, or both processes.Nevertheless, it is commonly accepted that dispersal limitation is a neutral stochastic process (Hubbell, 2001;Zhou and Ning, 2017).

Stochastic processes predominate the assembly of microeukaryotic communities
The neutral community model revealed the significant influence of stochastic processes on the assembly of microeukaryotic communities (Figure 10).Based on the robust assumption of ecological equivalence (Hubbell, 2005), this theory allows quantification of elusive processes-birth or death, dispersal, speciation, and ecological drift.The model has been widely applied to elucidate various ecological phenomena due to its simplicity and predictability (Heys et al., 2020;Jiao et al., 2020).The NCM explained a significant portion of microeukaryotic community variation across seasons, with over 63% of variation for each season and over 55% in total (Figure 10).This implies a stronger impact of stochastic processes on microeukaryotic communities compared to deterministic processes.The dominance of stochastic processes was reinforced by the results of VPA (Figure 9), which demonstrated that spatial factors contribute more to shaping microeukaryotic communities than environmental factors.These findings are consistent with previous research, for instance, Chen et al. (2019) revealed that microeukaryotic communities in river ecosystem were mainly shaped by stochastic processes, with 88.5 to 89.9% of the community variation explained by NCM.Similarly, in a study of Zou K. et al. (2021) investigating protist communities assembly in the Pearl River Delta with serious anthropogenic disturbance, NCM accounted for 51.24% ~ 75.82% of community variation.These findings collectively demonstrated that the assembly of freshwater microbial eukaryotic plankton was primarily driven by stochastic processes.
Variance partitioning analysis (VPA) is a valuable statistical approach extensively employed for evaluating the importance of environmental and spatial factors in the variation of microbial community ecology (Wu et al., 2020;Jing et al., 2023).Regrettably, our VPA models revealed that over 50% of the community variation remained unaccounted for (Figure 9).This observation aligns with several previous studies on microeukaryotic biogeography (Chen et al., 2019;Zhu et al., 2021), suggesting a limited impact of environmental and spatial factors on the microeukaryotic community.The large unexplained proportion could be due to influential factors not detected and included in the VPA.Furthermore, the inherent complexity of microbial ecological systems with intricate interactions, such as synergistic or antagonistic effects, poses challenges that may not be easily addressed by VPA (Meidute et al., 2008).It is also noteworthy that VPA tends to exaggerate the impact of spatial factors (Gilbert and Bennett, 2010).Consequently, enhancing the explanatory capacity of VPA may necessitate a broader scope of data acquisition, and the exploration of alternative modeling strategies, such as incorporating neutral community model, to ameliorate the inherent limitations of VPA.

FIGURE 1
FIGURE 1Seasonal and vertical variations in species richness (A) and Shannon diversity (C) of microeukaryotic plankton, and both fitted significant hump-shaped relationships with concentrations of chlorophyll a (B,D).
(2)  Promoting oxygenation to deeper hypoxic zones, stimulating the planktonic metabolic activities; (3) Facilitating the movement of microorganisms, increasing gene flow, genetic diversity, and the potential colonization of new habitats; and (4) Temporal changes in microbial community structure caused by the disruption of stratification, understanding these dynamics can help us gain insights into the resilience and adaptability of microbial communities.

FIGURE 8
FIGURE 8Drivers of planktonic microeukaryotic functional composition variation in Qianxia Lake.Random forest algorithm was used to determine the importance of environmental factors in explaining the variations of functional groups (A).The relative proportions of phototrophs, mixotrophs, and heterotrophs covaried well with several environmental gradients, including N:P ratio (B), Chl-a (C), and dissolved oxygen (D).Statistical significance is indicated as follows: *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 9
FIGURE 9Redundancy analysis (RDA) and variation partitioning analysis (VPA) depict the effects of environmental and spatial (PCNM) factors on microeukaryotic communities variations during spring (A,C) and autumn (B,D) in Qianxia Lake."E" represents relative contribution of environmental factors to community variation; "S" indicates relative contribution of spatial factors to community variation; "E&S" represents the shared explained variation; "Residuals" denote unexplained community variation.Statistical significance was assessed through permutational ANOVA (*p < 0.05 and ***p < 0.001).

FIGURE 10
FIGURE 10Predicted occurrence frequencies of amplicon sequence variants (ASVs) plotted against mean relative abundance, representing microeukaryotic communities in spring (A), autumn (B), and both seasons (C), respectively.Variations in proportions of neutral and non-neutral partitions in relation to microeukaryotic richness and abundance (D).The Nm value is determined by multiplying the metacommunity size (N) with immigration rate (m).The R 2 value indicates the fit adequacy to the neutral model.

TABLE 1
ADONIS analysis to assessed microeukaryotic communities variations using Bray-Curtis and Sørensen distances across seasons, and among depths and levels of chlorophyll-a.Significant differences in pair-wise comparison were bold.10.3389/fmicb.2024.1351772Frontiers in Microbiology 07 frontiersin.org