The Daqing River system (38°10′–40°102′ N, 113°39′′–117°34′′ E) is located in the central part of the Haihe Basin and crosses 64 counties and districts in the Shanxi, Hebei, Beijing, and Tianjin provinces. Its upper reaches are divided into two branches, the south and the north, which merge into Baiyangdian Lake (BYDL). The BYDL (38°3′–40°4′ N, 113°39′′–116°12′′ E), located in the northern part of the North China Plain, is a major part of the Daqing River water system (Zhao et al., 2011) ( Figure 1 ). The BYDL is a flood-relieving, stagnant, and comprehensively utilized lake connecting the Daqing River’s mountain areas and the plains. The Daqing River’s lower reaches enter the sea through the mainstream of Hai River in the north and Duliujian River in the south. In this study, we focused on the spatial distribution of macroinvertebrate communities in river-lake systems. The sample points were selected to cover the upstream, BYDL, and downstream regions, thereby reflecting the observed differences in the overall characteristics of the macroinvertebrate communities within the basin. A total of 152 samples were collected, 83 in 2019 and 69 in 2020, of which 69 were overlapping sample sites and 14 were not collected due to unforeseen circumstances.

A combination of type D trawls and modified Peterson mud harvesters were used to collect macroinvertebrate samples, and three replicate samples were collected at each sampling site (Zhang et al., 2021). The three replicate samples were mixed into one sample. After being sieved through a 60-mesh screen (250 μm) and removing any impurities, the samples were preserved with 75% alcohol and fixed with formaldehyde. A small amount of soda or borax was added to neutralize the formaldehyde to avoid corrosion of the calcareous shells, and insects and crustaceans were fixed with 75% alcohol in the sample bottles. The macroinvertebrates were identified to the lowest possible taxonomic level using appropriate identification guides under a microscope (Morse et al., 1994; Costas et al., 2018).

Physicochemical parameters, including water temperature (WT; ℃), dissolved oxygen (DO; mg/L), electrical conductivity (EC;μs/m), pH, and turbidity were measured using a YSI multiparameter probe (Yellow Spring, Ohio, USA). Water samples were collected from each sampling point using 500 mL plastic bottles and transported to the laboratory for water quality analysis. The ammonium nitrogen (NH₄-N), total phosphorus (TP), and total nitrogen (TN) contents were determined by UV spectrophotometry. The chemical oxygen demand (COD) was determined by potassium dichromate colorimetry. The measurements were conducted according to the relevant standards for water chemistry measurement in China (Huang et al., 1999). The environmental factors for the 152 samples under investigation are presented in Supplementary Table S1 . The water samples were collected before the macroinvertebrate samples to avoid any disturbance caused by the collection of the latter.

Spatial factors were employed as proxies for dispersal processes (Borcard et al., 2004). To create spatial factors, a Moran feature vector map (MEM) based on the Euclidean distances between sites was created according to geographic coordinates (x and y) (Blanchet et al., 2008). The principal coordinates of neighbor matrices (PCNM) is a method for the decomposition of spatial relationships among sampling sites, based on the diagonalization of a spatial weighted matrix (Bellier et al., 2007). For PCNM, the Euclidean distance matrix between the sampling points was calculated, and then the scale of the distance matrix was reduced. During this process, a Euclidean distance threshold needs to be set. Distances smaller than the threshold were retained, while distances larger than this threshold were set to four times the threshold. Quantitative ecologists have improved the traditional PCNM method, which can model any type of spatial structure. Thus, any similarity matrix can replace the Euclidean distance matrix to analyze the spatial structure, which is called distance-based Moran’s eigenvector maps (dbMEM). MEMs (i.e., spatial eigenvalues) can be used to model the spatial structure of biological assemblages at multiple scales, and they can represent dispersal processes, biogeographic effects, and environmental autocorrelation (Dray et al., 2012). In this study, the MEM vectors that were associated with positive eigenvalues and a significant Moran’s I were retained as they represented a positive spatial autocorrelation. The MEM was constructed using the “PCNM” package in R (Legendre et al., 2015).

All of the analyses were performed in R (version 4.3.2) unless otherwise indicated (R Core Team, 2022). The α-diversity index for each sample (i.e., the Shannon-Winner index) and β-diversity index (i.e., the “Bray-Curtis” distance) between the pairwise sampling sites were calculated using the “vegan” package. A principle coordinates analysis (PCoA) based on the “Bray-Curtis” distance was performed using the “vegan” package (version 2.6-4) to show the profile of the macroinvertebrate communities. Non-parametric statistical methods based on Adonis (Anderson, 2017) were employed to analyze the differences in macroinvertebrate community composition within rivers and lakes using the “vegan” package (Oksanen et al., 2019).

To address whether the DDR patterns of macroinvertebrate communities are variable with scale, the relationship between macroinvertebrate community similarity and geographic distance at the continuum river-lake system was investigated in 2019 and 2020. The study area was divided into two waterbody types (river and lake) and the distance decay relationship was calculated for each sampling point within the two spaces. The relationship between macroinvertebrate community similarity and geographic distance was calculated by linear least squares regression using the “vegan” package (Martiny et al., 2011).

To test the overall correlation between the two matrices of macroinvertebrate community and environmental factors and identify the potential driving factors of the formation of ecological patterns at the macro level, a mantel test was applied. However, a mantel test is highly dependent on the randomness and independence of the data. The importance of the effect of spatial and environmental variables on macroinvertebrate community variation was assessed separately based on redundancy analyses (RDA) and variance partitioning analysis. Variance Partitioning Analysis can determine the proportion of changes in community structure due to a given environmental factor. The contribution of each environmental factor and the interaction among the different environmental factors to the change in biomes can be calculated using sequencing analysis. However, the traditional varpart function in the “vegan” package can only handle four groups of explanatory variables. To improve normality and mean square deviation before the statistical analyses, the macroinvertebrate data were Hellinger transformed, and the environmental parameters were log (x+1) transformed except for the pH. For the RDA, spatial factors were calculated from the geographic distances using a PCoA of dbMEMs (Peres-Neto et al., 2006). The geographical factor (dbMEM) was determined using the “pcnm” function in the “vegan” package. To reduce multicollinearity between the spatial and environmental factors, variables with variance inflation factors less than 10 were retained in the analysis. The Euclidean ambient distance for each sample point in the two spaces was calculated to determine whether the environmental variables of the two water types were different.

To analyze the relative importance of the effect of the geographical distance and environmental factors on macroinvertebrate communities, a multiple regression on matrices (MRM) approach was employed. A MRM can consider multiple candidate models at the same time, evaluate each model, and select the model combination that can best explain the dynamic change in the vegetation community. However, the process of model selection is subjective. In this study, MRM was used for the different waterbody types (i.e., river and lake). The macroinvertebrate community dissimilarity matrix was calculated based on the abundance‐weighted UniFrac distance. A modified MRM code based on the R package “ecodist” (Goslee, 2007) was used to elucidate potential relationships between macroinvertebrate β-diversity and factors in the river and lake. The R2 values of the MRM model represent the total explanatory power of all the factors involved in the model, while the partial regression coefficient (i.e., b) represents the relative contribution of each factor. To remove covariant factors, non-significant factors were removed (Harrell, 2001). The MRM model was then performed again following the methods detailed in (Martiny et al., 2011).

To assess the importance of the effect of different processes on macroinvertebrate communities, the direct and indirect effects of spatial processes, physicochemical processes, and nutrients on macroinvertebrate communities were modeled using the “plspm” package (Liu et al., 2019). The corresponding influencing processes were represented by factors at different scales. The initial model included all possible pathways between the macroinvertebrate community composition, spatial geographic factors (MEM), physicochemical factors (WT, EC, DO, Turbidity, pH, and COD), and nutrients (TN, TP, and NH₄-N). A PCoA was used to summarize the variation in the distribution of the macroinvertebrate communities. The variables with loadings < 0.6 were removed and the goodness of fit index and R² were used to estimate the performance of the model. The “lavaan” package was used to perform partial least squares path model (PLS-PM) analysis (Mamet et al., 2019).

To ascertain the potential influence of stochastic processes on community assembly, the NCM (neutral community model) was used to predict the relationship between macroinvertebrate detection frequency and their relative abundance across the wider metacommunity (Sloan et al., 2006). The model presented here is an adaptation of the neutral theory adjusted to macroinvertebrate populations (Hubbell, 2001). In this model, Nm is an estimate of dispersal between communities, and it can be used to determine the correlation between occurrence frequency and regional relative abundance, with N representing the metacommunity size and m indicating the immigration rate. The parameter R² represents the overall fit to the neutral model (Sloan et al., 2006). All the fitting statistics were calculated with 1000 bootstrap replicates, and 95% confidence intervals were generated. The new mathematical framework of quantitative description of ecological stochastic processes that was proposed by Ning et al. (2019) was used, and the NST (normalized stochasticity ratio)() was employed to measure the importance of deterministic (homogenous and heterogenous) and stochastic (dispersal limitation and drifts) assembly processes. The assembly was considered to be dominated by deterministic processes (with a ratio of < 50%) or stochastic processes with a ratio of > 50% (Ning et al., 2019).

In both years, the river and lake macroinvertebrate communities comprised 72 taxa distributed across three phyla, six classes, and 15 orders ( Supplementary Table S1 ). Arthropoda was mainly distributed in the river, while Mollusca was widely distributed in the lake ( Supplementary Figure S1 ). Specifically, the relative abundance of Annelida and Arthropoda was lower in the lake than in the river, whereas Mollusca was higher ( Supplementary Figure S1 ). The macroinvertebrate α-diversity significantly differed between the river and lake in 2019. Specifically, compared with the lake, the observed Shannon and Simpson indexes were significantly lower in the river ( Supplementary Figure S1 ). Neither the Simpson Index nor the Shannon Index of the lake and river were significantly different in 2020. No significant difference in richness was observed inter-annually ( Supplementary Figure S2 ).

The PCoA results demonstrated a degree of dissimilarity between the samples, with the river and lake showing significant spatial variability, which was further supported by the non-parametric statistical approaches (R² = 0.09, P = 0.001, 2019; R² = 0.09, P = 0.001, 2020; Supplementary Figure S3 ). The PCoA analysis showed clear separation of communities by sampling location, with all the samples essentially grouped into two geographic clusters (i.e., river and lake; Supplementary Figure S3 ).

The community similarity of the lake was significantly higher than that of the river in 2019 and 2020 ( Figure 2 ). Significant DDRs were identified in both the river and lake macroinvertebrate communities, with the similarity of the macroinvertebrate communities significantly decreasing with increasing geographic distance based on Bray-Curtis similarity (P < 0.01, Figure 2 ). Moreover, between the two waterbody types, there were significant differences in the slopes of the DDRs. In the lake, the DDR slope was steeper, while the slope was flatter and the correlation was weaker in the river ( Figure 2 ). In addition, the DDR patterns of the macroinvertebrate community differed in the river between 2019 and 2020 ( Figure 2 ). Furthermore, significant DDRs were observed both upstream and downstream in the river ( Supplementary Figure S4 ).

Environmental distance decay was also significantly different in the different waterbody types in 2019 and 2020. In 2019, there was a significant environmental distance decay pattern for the macroinvertebrate communities in the lake, but the pattern was weaker in the river ( Supplementary Figure S5 ).

The Mantel test was employed to identify which environmental variables exhibited a significant correlation with the macroinvertebrate community. In 2019, pH showed a relatively weak relationship with the macroinvertebrate communities in the river (r <0.2, P < 0.01). Longitude (r > 0.2, P < 0.01) and TP (r > 0.2, P < 0.05) showed significant correlations with the macroinvertebrate composition in the lake. In 2020, WT had significant correlations with the macroinvertebrate community composition (r > 0.2, P < 0.05) in both the river and lake. Additionally, pH had strong and significant correlations with the macroinvertebrate community composition in the lake (r > 0.2, P < 0.05; Figure 3 ).

According to the VPA, the macroinvertebrate communities were influenced by both spatial and environmental factors. Macroinvertebrate communities in rivers were slightly more affected by spatial factors than environmental factors. Whereas, in 2020, macroinvertebrate communities in the lake were more affected by environmental factors than spatial factors, and, in 2019, they were more affected by spatial factors than environmental factors ( Figure 4 ).

In the MRM analyses, pH was the most important abiotic factor in explaining community variation in the river (b = -0.254, P < 0.05) in 2019. In the lake, the dominant factor driving changes in macroinvertebrate communities was geographical distance (b = -0.163, P < 0.05; Supplementary Table S1 ). In 2020, geographical distance and WT were the most important environmental factors in explaining community variation in the river (b = -0.314, P < 0.05; b = -0.522, P < 0.05). In the lake, the most important driver of changes in macroinvertebrate communities was WT (b = -0.488, P < 0.001; Supplementary Table S3 ). The varying importance of geographic distance and environmental factors reflected the different underlying variability across spatial-temporal scales in macroinvertebrate communities ( Supplementary Table S3 ).

We used PLS-PM analysis to explore the potential direct and indirect impacts of geographical factors, physicochemical factors, and nutrients on macroinvertebrate community composition (Supplementary Figure S6). The variation in the PCoA was explained by geographical factors, physicochemical factors, and nutrient variables. Spatial variables directly influenced macroinvertebrate community distribution, and they indirectly exerted influence via environmental factors. In 2019, the geographical variables, nutrients, and physicochemicals all had positive effects on the PCoA in the river, while they had negative effects in the lake. In 2020, the physicochemicals in the river directly affected the PCoA, while nutrients directly affected the PCoA in the lake.

The NCM accurately predicted the relationship between the occurrence frequency of the macroinvertebrates and variations in their relative abundance. The NCM explained 36.4%; 2019, river). 81%; 2019, lake). 69.6%;2020, river), and 83.1%;2020, lake) of community variance ( Figure 5 ). The NST results demonstrated that stochastic processes played a large role in shaping the macroinvertebrate community composition in the lake, but deterministic process dominated in the river (Figure 6).

The distance decay of similarity concept is inherent in many previous studies (Preston, 1962) and has been reported for various taxa and geographical settings (Whittaker, 1960; Soininen et al., 2007). In this study the waterbody type dependence of DDR patterns was investigated in a river-lake system. The results supported the hypothesis that the river zone decays at a slower rate when each spatial region of the Daqing River system was examined separately. The decay rates were more gradual in the river zone than in the lake. Furthermore, the results supported the other hypotheses to some extent, as the similarity of the macroinvertebrates in the lake decreased with distance. However, the similarity of the riverine macroinvertebrates did not differ significantly with increasing distance. The disparities among the different spatial ecoregions highlight the diverse factors that affect organisms at varying spatial levels (Maloney and Munguia, 2011). Some studies have shown that significant DDRs occur in macroinvertebrate communities (Astorga et al., 2012; Maloney and Munguia, 2011; Soininen et al., 2007). The relationship between macroinvertebrate community similarity and geographical distance may be stronger in lakes (smaller spatial scales). Our results are also consistent with the importance of geographic distance being limited to relatively small spatial scales and taxa with poor dispersal abilities (Astorga et al., 2012).

Our most noteworthy discovery was that the slope of the DDR became less steep as the spatial scale for the macroinvertebrate communities increased, which supports the hypothesis of Maloney and Munguia (2011). This suggested that macroinvertebrate community similarity decreased rapidly with the increase in geographic distance (particularly in the lake). This may be due to several ecological mechanisms. Lakes tend to have more homogeneous environmental conditions and lower habitat complexity than rivers, which allows for stronger effects of dispersal limitation and environmental filtering on community structure. In addition, the smaller the spatial scale of lakes, the greater the influence of local environmental factors, and the faster community similarity declines with distance. In addition, connectivity between lake habitats may be more homogeneous, leading to more consistent dispersal processes across the lake. These factors combine to influence distance decay patterns in lake macroinvertebrate communities. We also found that the slope of the DDRs was significantly more gradual in the river than in the lake. As macroinvertebrates are dispersed throughout the entire basin by water currents, the selective pressure exerted by the environment can significantly alter the community composition. Usually, the greater the disparity in the environmental conditions between localities, the more pronounced the variation in the macroinvertebrate communities (Heino et al., 2015b; Kärnä et al., 2015). This finding is consistent with other studies in phytoplankton, diatoms, and zooplankton (Bellier et al., 2014; Mazaris et al., 2010; Marquardt et al., 2018).

There may be several reasons for the waterbody type’s dependence on macroinvertebrates DDRs. First, the environmental heterogeneity was significantly larger at the river scale than at the lake scale ( Supplementary Figure S7 ). Second, the initial community similarity of adjacent sampling stations in the lake was high, which usually results in strong distance decay (Soininen et al., 2011). Third, the dispersal rate strongly regulates the compositional response of macroinvertebrate communities (Li et al., 2022; Zhao et al., 2024). The downstream areas were more impacted by human activities and exhibited greater homogenization compared with the headwaters. This led to a decrease in both taxonomic and functional diversity within the macroinvertebrate communities. Hence, human activities played a significant role in community shift, indirectly contributing to biogenic homogenization through alterations in physical and chemical properties, and subsequently influencing ecological processes (Ma et al., 2023).

Variance partitioning revealed substantial combined effects of the spatial and environmental variables. The Daqing River Basin spanned approximately 300 km. Therefore, the effects of spatial and environmental variables on the macroinvertebrate assemblages in the river were not surprising. The findings emphasize the significant impact of spatial variables on the variations in macroinvertebrate communities throughout the river ( Figure 4 ). In the river, macroinvertebrate communities with a greater geographical distance between them exhibited increasingly larger differences, which confirms the influence of spatial variables ( Figure 2 ). However, being surrounded or divided by farmlands, ditches, or paved roads related to agricultural activities may cause severe dispersal limitations for macroinvertebrates in the river (Guan et al., 2024).

A high degree of environmental heterogeneity was found within the lake and in the river, and both spatial and environmental factors influenced the community composition. The results of the VPA and MRM demonstrated that spatial factors exerted a slightly more pronounced influence on the structure of the macroinvertebrate communities than environmental factors. This is consistent with the findings observed in the Hun-Tai watershed (Zhou et al., 2020). The results of the VPA analysis in this study indicated that spatial factors played a significant role in shaping the macroinvertebrate communities in the river in both 2019 and 2020. In contrast, the lake region exhibited disparate outcomes between the two years, with environmental factors serving as the primary driver of community differences in 2020 (Spasojevic and Suding, 2012). This may be attributed to the high environmental heterogeneity observed in the lake in 2020. Geographical factors, physicochemical factors, and nutrients affected the macroinvertebrate community spatial distribution as found in other aquatic communities (Dong et al., 2024; Zhao et al., 2024).

The geographical factors (PCNMs) indirectly affected the distribution of the macroinvertebrate communities by altering the environmental factors. The macroinvertebrate communities in the river were not nutrition-restricted, while the lake had lower nutrient levels than the river ( Supplementary Figure S7 ). This may be due to the dilution effect of the lake on nutrients. Nutrient restriction affects primary productivity and, therefore, macroinvertebrate communities. To illustrate the effects of the spatial and environmental factors on macroinvertebrate communities’ spatial distribution, field observation should be conducted more frequently over a longer duration.

The NCM model exhibited a good fit for the majority of species, thereby substantiating the significant role of stochastic processes (i.e., dispersal and ecological drift) in the assembly of macroinvertebrate communities. The distance decay pattern of benthic macroinvertebrate communities in different water body types further confirmed the significance of stochastic processes ( Figure 2 ; Supplementary Figure S4 ). According to Hubbell’s Neutral Theory (Hubbell, 2001), due to dispersal limitation, community similarity is expected to decline with distance (Chase and Myers, 2011). Several studies have supported these findings. For instance, (Li et al., 2022) demonstrated that the benthic community structure in South Chinese bays was primarily driven by stochastic processes (R² = 0.90; NCM).

In this study, the R² value of the NCM parameter was slightly higher in the lake than in the river, and based on the Nm values, the macroinvertebrate dispersal between the sampling sites may have been higher in the lake than in the river (assuming a similar metacommunity size, Figure 5 ). This suggests that the effect of stochastic processes on the macroinvertebrates was stronger in the lake (Swenson and Enquist, 2009). Widespread dispersal is more likely to occur in highly connected water bodies (Tolonen et al., 2017), such as lakes (Cai et al., 2019) and main rivers (Brown and Swan, 2010). These results may be attributed to the higher habitat homogeneity in the lake compared with the river. High diffusion rates can reduce the effects of environmental selection and ecological drift to some extent. The wide distribution of the macroinvertebrates had a relatively small spatial impact on the community, suggesting that dispersal limitations had a weak effect on the spatial turnover of the macroinvertebrate communities in the lake.

Nevertheless, the NCM did not explain all the variation in the macroinvertebrate community, indicating that there might be other community assembly mechanisms, resulting in a non-neutral distribution. These other mechanisms include environmental selection and species interactions (Heino et al., 2015a). Stochastic processes establish the relationship between species occurrence and mean relative abundance (conforming to the neutral distribution), while environmental selection may modify this relationship. Species that are more adapted to local conditions have a higher frequency of occurrence irrespective of their relative abundance. This results in a non-neutral distribution of species occurrence. Another potential explanation is the varying migration or dispersal capabilities among species (Li et al., 2022). A higher occurrence frequency might indicate a higher dispersal ability than expected, and a lower occurrence frequency might indicate a lower dispersal ability than anticipated.

The relative importance of stochasticity in macroinvertebrate community building was further quantified by the NST based on a neutral model. The NST allows for more precise categorization of the relative contributions of stochastic and deterministic processes in benthic community construction and provides a more comprehensive perspective for understanding community construction mechanisms. The results suggested that deterministic processes may have a high influence in river macroinvertebrate community construction (Zhao et al., 2024), whereas construction may be more random in lake macroinvertebrate communities (Guan et al., 2024).

In summary, both ecological niche-related processes (environmental conditions and species interactions) and neutral processes (dispersal) were important for benthic community formation in this study, and both should be considered when exploring community assembly processes.

In this study, the diversity trends reflective of the river-lake system in temperate China. Additionally, the analysis was based on a two-year data set and does not include the full range of conditions that can occur in rivers and lakes, including wet periods, dry periods, and stable periods. Therefore, long-term data should be collected.

This study compared different regions and different years. Future research could use remote sensing techniques, such as satellite or unmanned aerial vehicle imagery. This information could provide large-scale, consistent, and repeatable measurements of land cover and water attributes, allowing for the validation of these research findings and determining the extent to which the observed patterns are generalizable. Remote sensing could be used to incorporate high-resolution hydrological conditions, such as wet and dry conditions, and it has reproducible frameworks (Valerio et al., 2024; Yue et al., 2023), which could enhance the understanding of the impacts of seasonal flow patterns on community composition. It could also be helpful for studying the dispersion of macroinvertebrates (Dong et al., 2021), which is strongly influenced by habitat landscape connectivity.