Geographic pattern of phytoplankton community and their drivers in lakes of middle and lower reaches of Yangtze River floodplain, China

Disentangling the relative contributions of deterministic and stochastic processes was critical to compressive understanding of underlying mechanism governing geographic pattern and assembly of phytoplankton community, while it was seldom performed in connected lakes under human pressure. Here, we investigated phytoplankton community pattern in relation to environmental and spatial factors over 81 lakes located in the middle and lower reaches of Yangtze River (MLYR) floodplain, where many lakes suffered from eutrophication and cyanobacterial blooms. A majority of MLYR lakes had higher phytoplankton abundance surpassing 107 cells/L and were dominated by common bloom-forming cyanobacterial genera, including Pseudanabaena, Microcystis, Merismopedia, Dolichospermum, Limnothrix, and Raphidiopsis. Phytoplankton community exhibited a striking geographical pattern both for taxonomic and functional compositions, while functional groups were less sensitive, and dissimilarity in communities displayed no significant increases with increasing geographical distance. Further, species richness explained much higher percentage of community variations than species turnover, indicating a reduced effect of environmental filtering of phytoplankton species with tolerance to similar environments in connected MLYR lakes. Both deterministic and stochastic processes governed assembly and biogeographic of phytoplankton community. Variation partition analysis showed that spatial factors exhibited greater influence on phytoplankton community compared to environmental variables. The stronger influence of spatial factors was further demonstrated by Mantel test and neutral community model. These findings indicate that deterministic and stochastic processes exhibited similar biogeographic patterns for phytoplankton community in MLYR lakes, but stochastic process was overwhelmingly dominated. Moreover, a large proportion of unexplained variation implies that complex interactions exist to shape assembly mechanism of phytoplankton community in MLYR lakes.


Introduction
Phytoplankton are a group of photosynthetic organisms with the size range from less than 2-200 μm, contributing with half of primary production on the earth and playing a pivotal Responsible Editor: Boqiang Qin Xiaochuang Li is responsible for managing all communication between the Journal and all co-authors, before and after publication. role in nutrient cycling, energy flow, and food web dynamics (Cardinale et al. 2002;Field et al. 1998). Phytoplankton are sensitive to environmental changes and have been widely demonstrated as an ideal bio-indicator for water quality assessment in a vast number of water bodies (Salmaso et al. 2006;Wu et al. 2019). In recent decades, anthropogenic inputs of excessive nutrients into lakes and reservoirs combined with global warming have resulted in recurrent blooms and further a rapid loss of biodiversity, ecosystem integrity, and stability in freshwaters (Huisman et al. 2018;Sala et al. 2000;Vörösmarty et al. 2010). Understanding the ecological process and mechanism driving phytoplankton diversity and biogeographic pattern is hence of great importance to provide baseline information to develop sustainable strategies for the management and conservation of biodiversity. A metacommunity is a set of local communities that are connected by dispersal of multiple potentially interacting species (Leibold et al. 2004). Traditionally, local environments (deterministic process) have been concentrated as the main drivers for phytoplankton community assembly from the aspects of bottom-up control by resources (light and nutrients) and top-down control by zooplankton and fish (Reynolds 2006). This could be ascribed that phytoplankton have small size and large population abundance that reduce extinction rate and increase dispersal probability and are assumed to be not limited by dispersal and recognized as cosmopolitans (Padisák et al. 2016). Contrarily, the neutral theory challenged this view by assuming that the interacting species are equivalent, and population dynamics are driving by random variation in births and deaths, of which dispersal ability is restricted by increasing geographical distance (stochastic process) (Hubbell 2001). A 25 years' survey in comparison with previous records of phytoplankton diversity covering 161 water bodies across 110,994 km 2 in Bulgaria found that 318 in 1393 infrageneric taxa exhibited restricted spatial distribution (Stoyneva 2016), while according to metacommunity theory, community assembly process can follow the neutral model or be explained by differences among species that related to their adaptation to local environmental factors and to their differential dispersal abilities (Rojo et al. 2016), suggesting that environmental and spatial factors act together on shaping community structure and diversity.
The metacommunity concept advances our understanding of the underlying mechanism that governs community assembly and helps to explain how organisms are dispersed in a given region through quantifying the relative contributions of environmental and spatial factors (Mihaljevic 2012), based on different theoretical paradigms (species sorting, patch dynamics, mass effect, and neutral model). Patch dynamics and mass effects are special cases of species sorting, emphasizing local processes, while neutral model assumes that traits shared among species show no difference and their population dynamics are driven by random losses and gains of species, focusing on regional processes (Winegardner et al. 2012). These four paradigms occupied not all of the inference of metacommunity theory for the mechanism structuring community assembly; meanwhile, they were not mutually exclusive (Brown et al. 2017;Winegardner et al. 2012).
In recent years, a growing number of studies have included space in the analysis but came to inconsistent conclusions regarding the relative importance of deterministic and stochastic processes in shaping geographical distribution pattern of phytoplankton communities. Some authors found that either spatial or environmental factors were the main drivers (Meier and Soininen 2014;Vanormelingen et al. 2008), while others demonstrated that both of these two factors played crucial roles in structuring community (Soininen et al. 2007). Nonetheless, some authors reported that neither environmental nor spatial factors explained a significant proportion of variations in phytoplankton community (Nabout et al. 2009). It was argued that inconsistent conclusions were reflected the differences in the spatial scale, the level of eutrophication, the strength of local environmental gradients, and the degree of connectivity among the studied freshwaters (Özkan et al. 2013;Soininen et al. 2011;Soininen 2014;Xiao et al. 2016). Taken the connectivity as an example, phytoplankton species were easily transported between connected waters, resulting in a homogenous community composition (Cottenie 2005). The connectivity reduced the filtering effects of local environmental conditions, while in isolated waters, dispersal was blocked, and the effects of local processes overrode that of regional processes in determining phytoplankton community structure (Padisák et al. 2010).
The middle and lower reaches of the Yangtze River basin are featured as the typical floodplain in which many shallow lakes connect with each other through intricate river networks (Yin et al. 2007). Particularly in the wet season, the frequent change of regional water level induces a mutual hydrodynamic exchange among these lakes, finally resulting in spatial destructuring and limnological homogenization (Bai et al. 2020;Thomaz et al. 2007). The floodplain lakes provide a good system to test the metacommunity theory for these lakes that are highly dynamic ecosystems with a wide diversity of habitats, as well as in their degree of connectivity (Ward and Tockner 2001). Moreover, a number of lakes suffered eutrophication and cyanobacterial blooms in MLYR (Chen et al. 2020); therefore, understanding the main drivers in structuring phytoplankton community in lakes under human pressure is vital for biodiversity conservation and sustainable watershed management.
Until now, rare studies have investigated geographical distribution pattern of phytoplankton community structure in the largescale MLYR floodplain lakes, especially considering the relative roles of both spatial and environmental factors. In the study, we disentangled and quantified the contributions of deterministic and stochastic processes in structuring phytoplankton community in 81 lakes of MLYR floodplain during flooding period. Our main objectives were to (i) uncover the geographical distribution and variation of phytoplankton community within and among MLYR floodplain lakes, (ii) identify the main drivers for phytoplankton community succession, and (iii) disentangle the quantitative contributions of deterministic and stochastic processes for variations in phytoplankton community. Taken together, this research will enhance our understanding on how phytoplankton communities were structured for the lakes with connectivity and human pressure and will provide baseline information for biodiversity conservation and sustainable watershed management in MLYR floodplain lakes.

Study area, sampling, and environmental and spatial variables
The middle and lower reaches of Yangtze River (MLYR) floodplain located in the subtropical climatic zone, with more than 600 lakes larger than 1 km 2 , distributed within the basin (Xu et al. 2017). The lakes played the pivotal role in both regional economic development and ecological stability (Cui et al. 2013). Heavy nutrients and other contaminant input from local and upstream areas have resulted in the deteriorated water quality, and most of these lakes were eutrophic and even ecologically unbalanced (Qin 2002).
Environmental and spatial variables, including physical (temperature), eutrophication indicators (SD, TN, and TP), lake morphological characteristics (depth, area, capacity), hydrodynamics (water retention time), the climate (precipitation), and spatial factors (PCNMs), were used to quantify the deterministic and stochastic selections in structuring phytoplankton communities in lakes along the MLYR floodplain. Surface water (0.5 m) was collected from 81 lakes within MLYR floodplain during August 2013 (Fig. 1). Between 3 and 30 sampling sites for each lake were designed depending on the lake area (Table S1). The eutrophication indicators (SD, TN, and TP) and water surface temperature (T) were averaged with monthly monitoring data on August 2013 from the China National Environmental Monitoring Centre and provincial environmental monitoring stations (Fig. S1, Table S1). The information on the lake morphological characteristics (depth, area, capacity) and the climate dataset (precipitation) were obtained from published research articles (Table S1). The hydrodynamics (hydraulic retention time, HRT) was defined: where capacity = lake area * lake mean depth and annual runoff = surface collecting area * precipitation * runoff coefficient (Shi et al. 2018). For the spatial component, the principal coordinates of neighbor matrices (PCNMs) analysis was performed with the vegan package to calculate a set of spatial variables based on the longitude and latitude coordinates of center of each lake (Oksanen et al. 2013).

Morphological identification and functional groups
For phytoplankton diversity analyses, 1 l freshwater was immediately fixed with 1.5% Lugol's iodine solution in situ and then transported to the laboratory for species identification (Izaguirre et al. 2016). Phytoplankton species were identified by 400 × magnification under microscope (Nikon, Japan), and at least 500 cells were counted for each sample. Morphological identification followed the nomenclature by Hu and Wei (2006). The species identified by morphology were further merged into functional groups using the criteria of Reynolds et al. (2002) and Padisák et al. (2009), which were grouped by similar ecophysiological adaptations to environmental constraints.

Community diversity and structure
Sampling sites with observed phytoplankton species < 10 were abandoned. Then, phytoplankton species abundance averaged from sampling sites within a lake obtained phytoplankton community composition and structure of each lake. Community diversity, including the number of species, Shannon-Wiener and Simpson diversities, and Pielou evenness were calculated using PRIMER v.7.0 (Clarke and Gorley 2015). Beta-diversity was estimated based on Bray-Curtis dissimilarity of community diversity and abundance. Beta-diversity was further partitioned into two components as species richness and species turnover that represented the balanced variation and abundance gradient, respectively. The partition was implemented using the bray.part function of the betapart R package (Jiao et al. 2017).
Non-metric multidimensional scaling (NMDS) analysis was used to reveal geographical pattern of phytoplankton community among lakes located in different provinces. Analysis of similarity (ANOSIM) was used to investigate significant differences in phytoplankton communities between set groups. The global R in ANOSIM analysis ranges from 0 to 1 that represents community separation degree between groups, with R = 0 indicating no separation and R = 1 as complete separation. NMDS and ANOSIM analyses were all performed with PRIMER v.7.0 (Clarke and Gorley 2015).

Relationships between phytoplankton community and environmental and spatial variables
To explore the potential impacts of environmental and spatial variables, the relationships between Bray-Curtis dissimilarity of phytoplankton community and Euclidean distance of local environments and geographic distance were analyzed based on Spearman's rank correlations. Before analysis, environmental variables and phytoplankton abundance were log(x + 1) and fourth root transformed, respectively, to improve homoscedasticity and normality for multivariate statistical analyses (Lepš and Šmilauer 2003).
The detrended correspondence analysis (DCA) was performed to decide whether redundancy analysis (RDA) or canonical correspondence analysis was used to investigate the relationships between phytoplankton community and environmental and spatial variables, depending on the longest gradient lengths. The longest gradient lengths were < 3 for both taxonomic compositions and functional groups of phytoplankton community, indicating that RDA was suitable for both communities. Before RDA analysis, environmental and spatial variables exhibited high variance inflation factor (VIF) ˃ 20 were eliminated to avoid collinearity among factors. A forward selection was performed to select those significant explanatory variables for further analysis (Blanchet et al. 2008). Then, variation partitioning analysis (VPA) was used to quantify the relative contributions of environmental and spatial variables in shaping community composition with adjusted R 2 coefficients based on RDA analysis . The explained contributions of environmental and spatial factors included environmental variables (E), pure environmental variables (E|S) that excluded contributions from spatial variables, spatial variables (S), pure spatial variables (S|E) that excluded contributions from environmental variables, and the combined effects of environmental and spatial variables (S ∩ E). The residual proportion represented unexplained variance by explanatory environmental and spatial variables.
Several studies based on simulation models have shown that VPA failed to correctly predict the environmental and spatial inference of community structure (Gilbert and Bennett 2010); thus, Mantel test and partial Mantel tests were conducted to assess the relationships between phytoplankton community variations and environmental/spatial variables. VPA and Mantel tests were combined to assess and verify the relative importance of environmental and spatial variables. Mantel tests were realized by the vegdist and mantel function in the vegan R package, and RDA and VPA analyses were performed with Canoco version 4.5.

Neutral community model
Sloan's community model was used to assess the potential importance of stochastic process for phytoplankton community (Sloan et al. 2006). The model predicts the relationship between the observed species frequency in a set of local communities and their regional relative abundance across the wider metacommunity and could reflect the adaptation of neutral theory suitable for large microbial populations (Hubbell 2001). In the model, the parameter R 2 determines the overall fit to the neutral model. The parameter Nm determines the correlation between occurrence frequency and regional relative abundance, with N and m describing the metacommunity size and immigration rate, respectively (Sloan et al. 2006). The neutral model was run in MicEco R package.

Taxonomic and functional diversity of phytoplankton community
A total of 397 infrageneric taxa from 8 phytoplankton phyla were identified in MLYR floodplain lakes, including Cyanophyta, Chlorophyta, Bacillariophyta, Euglenophyta, Cryptophyta, Dinophyta, Chrysophyta, and Xanthophyta. Chlorophyta was most diverse with 174 infrageneric taxa being detected, followed by Bacillariophyta (94), Cyanophyta (70), and Euglenophyta (27), with 32 infrageneric taxa in the other 5 phyla. Jiangsu Province exhibited the highest phytoplankton diversity, while Hubei and Jiangxi provinces presented the lowest diversity, basing on the evaluation of the number of species, Simpson and Shannon-Wiener diversity indices. As well, Hubei Province showed the lowest Pielou evenness index, whereas in Hunan and Anhui provinces, they displayed the highest evenness diversity (Fig. S2).
A contrasting variation of phytoplankton abundance was found in different lakes along the MLYR floodplain. In Jiangsu, Hubei, Hunan, and Jiangxi provinces, total phytoplankton abundance of a majority of lakes surpassed 10 7 cells/L or even higher than 10 8 cells/L; however, in Anhui Province, lakes exhibited a relatively lower abundance than 10 7 cells/L (Fig. S3). In Jiangsu, Hubei, Hunan, and Jiangxi provinces, Cyanophyta was overwhelmingly dominant over other phytoplankton phyla in all lakes, except some lakes in Hubei Province where phytoplankton abundance was lower than 10 7 cells/L, where Cyanophyta and Chlorophyta codominated in the lakes (Fig. S4). In Anhui Province, lakes were dominated either by Cyanophyta, Chlorophyta, Bacillariophyta, or combinations of these phytoplankton phyla (Fig. S4). In total, 44 infrageneric taxa showed wide distributions found in over 1/3 MLYR lakes. Among these, Pseudanabaena, Microcystis, and Merismopedia were the most dominant for they had a mean abundance of 10 7 cells/L over the whole 81 lakes, followed by Dolichospermum, Limnothrix, and Raphidiopsis that had a mean abundance of 10 6 cells/L (Fig. 2a).
Altogether, 29 functional groups (FGs) were identified in MLYR floodplain lakes. Similar to phytoplankton diversity of taxonomic compositions, Jiangsu Province exhibited the highest functional diversity, while Hubei and Jiangxi provinces presented the lowest diversity. For Simpson evenness index, we came to the same result that Hunan and Hubei provinces showed the highest and lowest evenness diversity, respectively (Fig. S5). A total of 19 FGs exhibited wide distributions in over 1/3 MLYR floodplain lakes; among these, FGs belonging to Lo, S1, and M showed the high mean of phytoplankton abundances over 10 7 cells/L, and FGs attributed to SN, J, H1, MP, T, D, Tc, P, and S2 had the mean abundances over 10 6 cells/L in MLYR floodplain lakes (Fig. 2b).

Geographical distribution pattern of phytoplankton community
Lakes within five provinces showed a more or less separation of phytoplankton communities, both for taxonomic and functional compositions (Fig. 3). The phytoplankton community in Hubei lakes distributed randomly and had a low separation from that in lakes of other provinces. This was demonstrated by ANOSIM analysis that community separation was significant compared with lakes from different provinces, except comparisons between Hubei and other provinces; moreover, this was particularly evident for functional groups (Table 1).
Dissimilarity in phytoplankton community increased with increasing geographical distance ( Fig. 4a-b), while the positive correlations were significant for taxonomic composition (R = 0.098, P = 0.024) but not for functional groups (R = 0.061, P = 0.114). Beta-diversity partitioning revealed that species richness (nestedness) rather than species replacement (turnover) contributed to the majority of the beta-diversity and drove the community changes over geographical distance, for both taxonomical and functional compositions (Fig. 4c-d). Species richness displayed a significant increase with distance for taxonomic composition (R = 0.098, P = 0.024), while their correlations were not evident for functional composition (R = 0.042, P = 0.193).

Environmental and spatial factors related to patterns of phytoplankton community
We explored the balance between deterministic and stochastic selections in structuring biogeographic pattern of phytoplankton community in MLYR floodplain lakes. Phytoplankton community displayed a significant distance decay pattern with increasing geographic distance (Fig. 4a), indicating that spatial factors played the important role in structuring phytoplankton communities. Compared to taxonomic compositions, the functional groups of phytoplankton community showed a less evident distance decay pattern (Fig. 4a-b).
Inferred from RDA ordination plot, only spatial factors (PCNMs 2-5) were significantly related to taxonomic compositions (Fig. 5a), while for functional compositions, both environmental (TP and precipitation) and spatial (PCNMs 3 and 5) contributed to the significant changes of communities (Fig. 5b). The Mantel tests further demonstrated that spatial effects exhibited more influential impacts than environmental variables in structuring phytoplankton community, both for taxonomic and functional compositions (Table 2). PCNM1 and PCNM2 showed significant correlations with phytoplankton communities, while all of environmental factors had minor influences. At provincial scales, local Mean relative abundance of taxonomic (a) and functional phytoplankton diversity (b) in MLYR floodplain lakes. In (a), only infrageneric taxa occurred in over 1/3 lakes are shown; taxa belonged to Cyanophyta, Chlorophyta, Bacillariophyta, Euglenophyta, and Cryptophyta are marked in blue, green, crimson, pink, and gray colors.
Each of the squares at the bottom represents one lake, and the colors are corresponding to the disparity in total phytoplankton abundance. Numbers in the right indicate the occurrences of taxa shown in the 81 studied lakes environmental variables played the pivotal role in structuring phytoplankton communities especially for those in Jiangsu and Hunan lakes, whereas spatial selections displayed less influential impacts (Table 2).
VPA analysis revealed the similar results for taxonomic and functional compositions of phytoplankton community changes explained by environmental and spatial factors. For taxonomic compositions, the pure spatial variables (24.0%) explained a much higher proportion of community variation than that by pure environmental variables (3.0%) (Fig. 6a). As well, the pure spatial variables (23.5%) exhibited a higher explanation than by pure environmental variables (3.1%) for functional composition (Fig. 6b). The shared explanations by jointly environmental and spatial factors were equal to purely environmental variables. The strong effects of spatial factors were further demonstrated by the neutral community model (NCM). NCM assessed the potential importance of neutral processes of phytoplankton community and explained a large fraction of variations both for taxonomic (R 2 = 0.834) and functional (R 2 = 0.732) compositions (Fig. 7).

Phytoplankton community exhibited evident biogeographic patterns in MLYR floodplain lakes
With the rapid population increase and economic growth in MLYR floodplain, intensive anthropogenic activities brought excessive nutrients into lakes that resulted in water quality deterioration (Feng et al. 2019;Luo et al. 2019). Using MERIS and OLCI remote sensing images to retrieve chlorophyll-a concentrations, all 50 examined large lakes in MLYR floodplain presented high probabilities of eutrophication (Guan et al. 2020). This was further demonstrated in our study that a majority of MLYR floodplain lakes exhibited a much higher phytoplankton abundance surpassing 10 7 cells/L or even higher than 10 8 cells/L (Fig. S3). In these lakes, the common bloom forming cyanobacteria Pseudanabaena, Microcystis, Merismopedia, Dolichospermum, Limnothrix, and Raphidiopsis were predominant (Fig. 2a), suggesting that MLYR shallow lakes would be highly susceptible to cyanobacterial blooms. These dominant genera were affiliated to FGs as Lo, S1, and M that were favored by turbid, eutrophic, and high irradiance conditions (Fig. 2b), which were characteristics of shallow lakes (Reynolds et al. 2002).  The MLYR floodplain constitutes a network of shallow lakes with different connectivity degrees and allowed the exchange of hydrodynamics via watercourses in the wet season, resulting in an increased habitat similarity (Devercelli et al. 2016). Similarly, local environments exhibited similarities among lakes from different provinces in MLYR floodplain (Fig. S6). Different species had their own properties, such as physiological tolerances and dispersal abilities, and the similar environments might result in a final homogenization of phytoplankton communities. However, we found a striking geographical pattern of phytoplankton community among lakes of different provinces in MLYR region (Fig. 3). Though flooding generated the similar environments in connected lakes, yet it took time to facilitate dispersal among lakes; thus, spatial patterns existed in the previous period were still present during the flooding period (Rojo et al. 2016). Furthermore, Euclidean distance of environments showed significant increases with geographical distance (Fig. S7). The MLYR lakes were connected to some degrees, while inland lakes were separated by dry land that led to hydrodynamic gradients in these lakes (Bai et al. 2020). Additionally, a large number of dams built for reservoirs and flood control also acted as the barriers for the exchange of hydrodynamics among lakes. By the year of 2009, over 47,000 dams have been built in the Yangtze River basin . Contrarily, bacterioplankton community tended to cluster together in the lakes with small distance in MLYR floodplain (Bai et al. 2020). Bacterioplankton in Jiangsu lakes shared high similarities but exhibited significant dissimilarities from those in Anhui, Jiangxi, and Hubei lakes. The discrepant results were explained that decreasing size increased the dispersal rate of passively dispersing species; thus, larger phytoplankton exhibited more spatial limitation than by smaller bacterioplankton (Beisner et al. 2006;Mazaris et al. 2010). These complex mechanisms resulted in geographical separations along MLYR connected lakes in different provinces.
Beta-diversity linked local and regional diversity by measuring the amount of species dissimilarity between communities and played a pivotal role in disentangling the Fig. 4 Spearman's rank correlations between Bray-Curtis dissimilarity of phytoplankton community and geographic distance, including total beta-diversity, and the turnover (red dots) and richness (green dots) components of beta-diversity for taxonomic (a, c) and functional compositions (b, d) of phytoplankton community in MLYR lakes various processes that shaped biodiversity patterns across scales (Viana et al. 2015). Beta-diversity could be further divided into two components: species richness represented that local communities were formed by subsets of species from regional communities, while species turnover indicated that one species was substituted by another or the changes of its relative abundance in the community among sites (Ulrich et al. 2009;Yang et al. 2018). In MLYR lakes, beta-diversity partitioning revealed that species richness rather than species turnover governed the assembly and biogeographic pattern of phytoplankton community structure (Fig. 4c-d).
As shown in isolated freshwaters and freshwaters with contrasting environmental gradients, spatially structured environmental filtering induced the high species replacement, mediated by species properties adapted to different environments (Vanormelingen et al. 2008;Yang et al. 2018). Among interconnected lakes, matters and hydrodynamics were exchanged through watercourses that resulted in increased habitat similarity (Thomaz et al. 2007); in addition, hydrodynamics exhibited not significant differences in the floodplain lakes (Bai et al. 2020) and reduced the effect of environmental filtering of phytoplankton species with the tolerance to these similar environments. Moreover, spatial isolation could increase the relative importance of species richness and reduce the contributions of species turnover (Bender et al. 2016). The studied 81 connected lakes spanned an intermediate scale over 1000 km, highlighting the importance of spatial factors in structuring phytoplankton communities in MLYR lakes, and their relative contributions were further discussed below.

Stochastic processes contributed over deterministic processes in governing phytoplankton community
Metacommunity theory suggested both environmental and spatial factors shape community structure and diversity (Leibold et al. 2004). Due to their small size and rapid and complex responses to environmental changes of phytoplankton species, it was difficult to divide the relative contributions of environmental and spatial variables (Etienne and Alonso 2007). The geographical distance has been used for measuring the dispersal ability in shaping microbial communities  (Melo et al. 2011), and the relative importance of dispersal could be further demonstrated by the distance decay pattern between community and geographical distance . Phytoplankton community displayed a significant distance decay pattern with increasing geographic distance (Fig. 4a), strengthening the important role of spatial factors in structuring phytoplankton communities in MLYR lakes. Compared to taxonomic compositions, the functional groups of phytoplankton community showed a less evident distance decay pattern (Fig. 4a-b). Functional groups were composed of set of species that shared common physiological and ecological traits (Padisák et al. 2009;Reynolds et al. 2002), which reflected more sensitive responses to environmental gradients than by taxonomic compositions, thus resulting in much higher dispersal limitations with increasing distance scales. It could be expected that functional groups exhibited a more evident distance decay pattern than by taxonomic compositions in MLYR lakes, while we got the adverse results ( Fig. 4a-b). Connected lakes in MLYR lakes exchanged matters and organisms that resulted in the increased habitat similarity, and a homogenization of phytoplankton community might complicate these processes.
As well, the neutral model has been widely used to quantify the relative importance of stochastic processes (Sloan et al. 2006) and could correctly explain microbial community assembly and pattern in diverse environments (Roguet et al. 2015). The neutral model successfully explained a large proportion of variations in phytoplankton community and further demonstrated the pivotal role of spatial factors in governing community assembly (Fig. 7). As we expected, the explained variations of taxonomic (R 2 = 0.834) compositions were higher than that by functional groups (R 2 = 0.732), further demonstrating that spatial variables had more influential effects on taxonomic compositions than functional groups.
Previous studies concluded that both deterministic and stochastic processes governed the assembly and biogeographic pattern of phytoplankton community in aquatic systems (Devercelli et al. 2016;Izaguirre et al. 2016;Rojo et al. 2016;Soininen et al. 2007;Yang et al. 2018). It was argued that the relative role of deterministic and stochastic processes was related to the sampling scale (Izaguirre et al. 2016). At small scales (few kilometers), phytoplankton with highly passive dispersal abilities were easily transferred from one site to another, and local environments acted as the only filters for species. With the spatial scales increasing, geographical distance as barriers blocked their dispersal and led to a combination of both environmental and spatial constraints. At a wide spatial scale (tens of thousands of kilometers), the effects of regional processes overrode that of local processes in governing phytoplankton community. The MLYR lakes spanned an intermediate scale over 1000 km, and bias was obtained by Mantel tests that only spatial factors influenced the geographic pattern of phytoplankton communities (Table 2). Consistent with previous assumptions (Izaguirre et al. 2016), at provincial scales (small scales), local environments predominated over spatial variables in structuring phytoplankton communities ( Table 2).
The biased results at regional scales in MLYR lakes were also observed in RDA analysis that only spatial variables were significantly related to taxonomic compositions of phytoplankton (Fig. 5). VPA analysis further demonstrated that spatial variables contributed to a larger proportions over local environments in structuring phytoplankton communities (Fig. 6). Several explanations could account for the overwhelming superiority of spatial variables. First, as mentioned above, the studied 81 connected lakes spanned an intermediate scale over 1000 km in the MLYR floodplain, indicating a combination of both stochastic and deterministic processes involved in the assembly of phytoplankton community structure (Soininen 2012). Second, shallow lakes in MLYR floodplain connected in the wet season and exchanged hydrodynamics and organisms through watercourse among these lakes. The connectivity reduced the filtering effects of local environmental conditions and finally resulted in a homogenization of phytoplankton communities in MLYR floodplain lakes (Cottenie 2005). Third, the disordered interconnectivity among lakes without hydrological bias increased the random distribution of species distribution (Bai et al. 2020), further supporting a wide distribution of phytoplankton species in MLYR floodplain lakes. Last but not least, the stochastic process overwhelmed deterministic process when abundant bacterial species were generalists in various habitats (Liu et al. 2015). In MLYR floodplain lakes, 44 infrageneric phytoplankton taxa showed wide distributions in over 1/3 lakes (Fig. 2a), suggesting that these generalists contribute to the overwhelming superiority of stochastic process over deterministic process in governing geographic pattern of phytoplankton community in MLYR floodplain lakes.

Conclusion
In the study, both deterministic and stochastic processes drove the geographical pattern of phytoplankton community structure in MLYR floodplain lakes, although their relative contributions were different. However, spatial variables explained more community variation than by environmental factors for connected lakes in MLYR floodplain exchanged matters and organisms and resulted in the increased habitat similarity and homogenization of phytoplankton community. Though important local environments (TN and TP) were included, the explained variation by environmental variables was much low. Therefore, to protect MLYR floodplain lakes under human pressure, both local environmental conditions and the effect of regional phytoplankton species pool on local community should be considered. Our findings enhance the mechanistic understanding of geographic pattern of phytoplankton community among interconnected lakes under human pressure and have important implications for biodiversity conservation and ecosystem function in these connected and eutrophicated lakes.