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 vacuum-filtered 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₃-N), nitrite (NO₂-N), ammonium (NH₄⁺-N), total nitrogen (TN), and soluble reactive phosphate (PO₄³⁻) 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).

We extracted DNA using the PowerSoil^® DNA Isolation Kit (MoBio, Carlsbad, CA), following the manufacturer’s protocol. The assessment of DNA concentration and quality was performed with an NC2000 spectrophotometer (Thermo Scientific, United States). PCR amplification targeted the V4 region of microeukaryotic 18S rRNA gene, employing the specific primers 565F (5′-CCAGCASCYGCGGTAATTCC-3′) and 948R (5′-ACTTTCGTTCTTGATYRA-3′) (Stoeck et al., 2010). PCR reactions followed this procedure: 98°C for 2 min, followed by 25 cycles of 98°C for 15 s, 55°C for 30 s, and 72°C for 30 s; and a final extension at 72°C for 5 min. Paired-end (2 × 250 bp) sequencing was conducted using the Illumina NovaSeq platform (Illumina, United States) at Personalbio company (Shanghai, China).

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² 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).

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₁₀(χd) = (−2z) * log₁₀(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.

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² database, were converted into ecological traits, while retaining their read proportions within a given community. In this study, the trophic traits were simplified and categorized into phototrophy, heterotrophy, and mixotrophy (see Supplementary Table S4). The ASV richness and abundance for these three functional groups were determined by summing the taxonomic groups within each corresponding functional category.

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² were assessed using the A3 package in R (Fortmann-Roe, 2015).

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² 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.

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).

The physicochemical variables determined are presented in Supplementary Table S1. Of the 12 environmental parameters, more than half demonstrated a distinct seasonality. The water temperature ranged from 19.1 to 27.3°C in spring and from 26.4 to 29.9°C in autumn (p < 0.001). Significant differences were observed in Chl-a concentration, which was higher in spring (mean ± SE, 5.901 ± 0.502 μg/L) than in autumn (1.895 ± 0.203 μg/L; p < 0.001). Similarly, levels of DO (9.598 ± 0.138 vs. 8.069 ± 0.078 mg/L), NO₂-N (0.004 ± 0 vs. 0.003 ± 0 mg/L; p < 0.001), pH (9.435 ± 0.133 vs. 8.694 ± 0.079; p < 0.001), and TOC (70.208 ± 2.576 vs. 58.833 ± 0.746 mg/L; p < 0.001) showed spring-autumn distinctions. However, concentrations of TP (0.421 ± 0.011 vs. 1.350 ± 0.027 mg/L; p < 0.001) and PO₄-P (0.124 ± 0.009 vs. 0.153 ± 0.010 mg/L; p < 0.001) were relatively lower in spring due to increased nutrient uptake by high abundance of phytoplankton. In terms of depths, environmental factors were relatively stable, although the concentrations of DO, pH, and temperature significantly decreased with depth, while TN tended to increase with depth (Supplementary Table S2).

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² ≥ 0.115, p ≤ 0.033), and the species richness appears to peak at the concentrations of Chl-a of 4.0 μgL⁻¹ (Figures 1B,D).

The ASV number were mainly affiliated with phylum Chlorophyta (on average 27%), followed by Dinoflagellata (14%), Chytridiomycota (11%), Ciliophora (9%), Ochrophyta (9%), Cercozoa (4%), Cryptophyta (4%). Seasonal variations were noted, with higher occurrences of ASVs in the spring compared to the autumn for Ciliophora (11% vs. 7%), Ochrophyta (11% vs. 7%), and Cercozoa (5% vs. 3%), but the opposite for Dinoflagellata (8% vs. 20%) (Figure 2A).

The reads of microeukaryotes exhibited dominance by Chlorophyta, accounting for a relative abundance of 37.5%, followed by Fungi at 36.3%, mainly composed of Chytridiomycota and Basidiomycota. Other notable phyla included Alveolata (14.2%), Stramenopiles (5.6%), and Cryptophyta (2.9%). Within Alveolata, Ciliophora comprised 6.6% and Dinoflagellata constituted 6.4%. Stramenopiles were predominantly represented by Ochrophyta (2.9%) and Bacillariophyta (2.3%). Reads proportions of other taxa (e.g., Haptista, Cercozoa, Amoebozoa, and Opisthokonta) averaged less than 1% (Figure 2B).

Additionally, the ASV richness in Dinoflagellata, Basidiomycota, Oomycota, and Perkinsozoa exhibited significant hump-shaped relationships with Chl-a, except for Chytridiomycota and Bacillariophyta, which formed negative relationships with Chl-a (0.083 ≤ R² ≤ 0.357, p ≤ 0.044; Supplementary Figure S2).

The NMDS plot and PERMANOVA analysis indicated a distinct separation in microeukaryotic community compositions between spring and autumn (R² = 0.302, p = 0.001). However, samples from different depths did not exhibit clear differentiation (R² = 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² ≥ 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² ≥ 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).

Comparing the proportions of major microeukaryotic taxa in reads between seasons revealed that Chlorophyta (41.8 vs. 30.8%), Haptista (1.03 vs. 0.27%), and Choanoflagellida (0.40 vs. 0.05%) were significantly more prevalent in autumn than in spring. In contrast, Basidiomycota (0.28 vs. 9.34%), Dinoflagellata (3.35 vs. 11.1%), Apicomplexa (0.10 vs. 0.63%), and Mesomycetozoa (0.08 vs. 0.27%) exhibited an opposite trend (Figures 4A–G; Supplementary Table S5). Regarding the depth, the variations were mainly attributed to Ascomycota, Cercozoa, and Cryptophyta (Figures 4H–J; Supplementary Table S5).

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⁻¹ for low level; 2–5 μg L⁻¹ for intermediate level; and 5–11 μg L⁻¹ 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).

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 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² = 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₂-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² = 0.097, p = 0.016), whereas heterotrophs exhibited a significant and negative correlation with N:P ratio (R² = 0.157, p = 0.002). Nevertheless, the correlation between mixotrophs and N:P ratio was weak and non-significant (R² = 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).

RDA demonstrated variations in the percentage of explained variance for microeukaryotic communities by the first two axes across seasons, 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₂-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 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 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.

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₂-N, PO₄-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₂-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); (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.

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⁻¹) 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.

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 decoupling between taxonomic and functional diversity contrasts with the findings of Ramond et al. (2019), who indicated a strong association between the functional diversity of marine protist communities and their taxonomic diversity across different size fractions. The discrepancy might be ascribed to their transcontinental survey (~3,800 km) in contrast to our localized one (40 km). Studies on microbial biogeography and large-scale environmental surveys have documented higher diversity with larger sampling scales, potentially leading to increased functional diversity (Martiny et al., 2006). Additionally, in comparison to the various traits annotated by Ramond et al. (2019), we acknowledge limitations in simplifying to three functional traits.

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₄-P (on average 0.14 mg/L) but low levels of DIN (dominated by NO₃-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).

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).

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.