Linking the Community and Metacommunity Perspectives: Biotic Relationships Are Key in Benthic Diatom Ecology

: The ecology of benthic diatoms is scarce in diatom reviews, and it seems that the loss of interest in their local ecology (populations–communities) coincides with an increase in metacom-munity studies. We include a review of the latter to highlight some unresolved issues. We aim to demonstrate the relevance of local population–community ecology for a better understanding of the metacommunity by addressing gaps such as the relevance of biotic relationships. We analyzed 132 assemblages of benthic diatoms from two neighboring catchments, with varying altitudes, lentic and lotic waters and substrates. Population–community features (e


Community and Metacommunity of Benthic Diatoms: The Scale
Our purpose is to demonstrate the relevance of local population and community ecology for a better understanding of the metacommunity. Related to this, we highlight the role of biotic relationships as controlling factors which seem to be neglected in metacommunity studies. Take, for example, the case of benthic diatoms (BDs). The emergence and development of BD metacommunity analysis seems to coincide with the preclusion of more locally focused ecological studies (i.e., population and community) on BD. We advocate for a recovery of local ecology approaches, which are necessary to understand that diversity is the outcome of assembly dynamics, as recognized for plankton [1,2]. Thus, we suggest that this approach could be the bridge to a better understanding of the structure and dynamics of metacommunities (Table 1).
The main objectives in the study of BD metacommunities (Table 1) are usually to describe features such as beta diversity (hereafter βD), the underlying explanatory models [30,31] and the variance explained by their possible controlling factors (i.e., variance partitioning techniques; [32]). These BD studies are, with few exceptions, carried out over a large area and with a very different grain and lag (in the sense of [33]; Table 1). It is surprising that BD studies use the same scales as those used in studies of other (larger) aquatic organisms, such as fish [6], and that kilometric distances are considered "fine Table 1. Main results reported for metacommunity analyses of benthic diatom (BD) communities worldwide. The lag or spacing (km) was estimated as the square root of geographical extent scaled by the number of sampled sites. ANOSIM: ANOVA of similarities; CCA: canonical correspondence analysis; DBMEM: distance-based Moran's eigenvector maps; DCA: detrended correspondence analysis; DBRDA: distance-based RDA; LCBD: local contribution to beta diversity; LM: lineal model; MRM: multiple regression model; NMDS: non-metric multidimensional scaling; PCNM: principal coordinates of neighbor matrices; RDA: redundancy analysis; SAD: species abundance distribution; SCBD: species contribution to beta diversity. RDA on spatial and environmental factors and taxa

Site
The environment plays the most important role in structuring stream BDA, but spatial factors also explain some variation in diatom distribution, especially at the more coarse scale (i.e., continental) [4] Mesta river (Bulgaria) 5.0 × 10 3 (6.5) Cobble (whole range) Taxa relative abundances RDAs on environmental, temporal and spatial factors and taxa All three independent matrices explain variability in BDA [5] River Viaur (France) While richness and beta diversity of streams are related to the regional environment, those of lakes are related to spatial measurements; differential hydrological connectivity is the key factor of these diatom variables [19]  Genetic diversity was higher than morphodiversity; alpha and beta diversity responded differently to lake depth [29] Therefore, when looking for explanations of the assemblies of BD metacommunities, we should also remember that environmental factors can control BD spatial distribution at very fine scales, for example, within a pond [35] or in a few square meters of a stream [3]. To address these questions regarding the analyzed space, we used intermediate scales compared to those mentioned in Table 1. In Bailey's terminology [36], these would be a set of sites analyzed at a local scale (microecosystems or local communities with samples separated by less than 2 km) included in a smaller landscape mosaic of an intermediate scale (with a total area of over 3000 km 2 ). This approach can be effective in disentangling mechanisms acting at different organizational levels of BD (populations, communities and metacommunities).

Drivers of BD Distribution: Environmental Factors, Habitat and Substrate
Environments that are stressful for different taxa (e.g., extreme values of factors, high disturbance) can act as a filter (i.e., only populations of some adapted species can grow in such conditions, and therefore beta diversity will be low). On the contrary, non-stressful environments attain a higher βD thanks to assembly processes due to the interacting species. The combination of species generates alternative communities and, therefore, a greater βD [2,37]. For years, there have been warnings that the idea of limiting factors in phytoplankton prevented a thorough interpretation of possible alternative mechanisms, such as biotic relationships (coexistence and exclusion) in the assembly process [2,38,39]. The relevance of biotic factors, such as the interaction between BD species, as a possible local control was reported decades ago [40,41], and more recently some authors have related them to community composition and spatial distribution [42,43]. However, these interactions have been neglected in BD distribution studies. The larger the area studied, the more likely it is to include environments with extreme values of abiotic factors, such as for conductivity or nutrient concentrations (Table 1), and the metacommunity will then be partly structured by this abiotic environment. In fact, the role of such factors (e.g., nutrients) can be concomitant with other factors, resulting in inconclusive patterns of response [4,23]. Including BD species abundances as potentially relevant environmental factors for other BD populations is a necessary step to elucidate if controlling factors are only abiotic or if they are also biotic. In addition, the effects of herbivores, such as fish and the allelopathic effect of macrophytes, on BD may also determine their community structure [35,[44][45][46]. Therefore, biotic factors should not be overlooked in metacommunity studies as they can exert additional local control in BD predictive models at different scales [10,47]. Disentangling the role of biotic relationships as drivers of populations, communities and metacommunities is an important goal of this study. To do so, we focus on an area where local abiotic conditions are neither extreme nor disturbed.
There are other interesting pieces to this puzzle about the causes of BD species distribution that can be addressed from a local-intermediate scale approach and that we analyze in this study: the effect of substrate, habitat and altitude. Many studies have relied on analyses of BD abundances growing on easily quantifiable substrates in terms of surface area, such as rocks and cobbles ( Table 1). The issue of substrate influence on organismic growth has been a long-standing topic in diatom ecology [40,[48][49][50], but it is far from being solved [51]. Very few studies address the comparison of the metacommunity of stagnant water BD with that of streams (i.e., [52,53]), despite the fact that the physical features of one or the other habitat may affect community connectivity, environmental homogenization or disturbance intensity [52,54]. Altitudinal effects on diatom communities have also been a preferred topic for diatom ecologists for many years [12,15,55], but the results are, as yet, pretty inconclusive [47].

Goals of This Study
Our main goals in this paper are to disentangle the differences in environmental factors controlling populations, communities and metacommunities of BD and to test the relevance of biotic relationships and the importance of their effects at the local scale. We will delve deeper into these questions using a BD field-based study to establish the main mechanisms and processes of assembly as a "bridge" to better explain the metacommunity structure. In addition, this will provide more information about the ecology of this group of organisms that has been neglected in previous reviews which have mainly focused on planktonic diatoms [40,[47][48][49]56,57].
Based on the issues mentioned here, we hypothesized that in a BD landscape of intermediate extension, which favors little variability in the intensity of disturbances (except for the obvious lotic and lentic differences): (i) α diversity will be independent of environmental factors and favored in undisturbed systems; (ii) relationships between BD populations will become very important for the assembly process, with exclusion (overdispersion) being the most relevant; (iii) βD will be higher if environmental stability favors alternative assemblages; and (iv) the metacommunity will be structured mostly by species sorting effects, particularly if we include the biotic factors as an environmental component.
The study was carried out in central Spain (Mediterranean semi-arid climate), thus widening the geographical scope of BD metacommunity knowledge, since the majority of such studies have, so far, been carried out in cold temperate environments (Table 1). We worked in two neighboring catchments with pristine ecosystems within an altitude range of 800 m.

Field Work
Sampling was carried out in a selected area (3588 km 2 ) in the highest, bordering part of two large catchments (the Júcar and Tajo rivers), a mountainous region of centraleastern Spain called Serranía de Cuenca ( Figure 1). The area is largely depopulated (3.5 inhabitants/km 2 mostly living in small towns and villages) due to migration to eastern Spanish cities since the 1960s [4]; these trends have not changed recently, and tourism is still scarce (more information on its geography in Supplementary Materials). Coordinates and more information on these locations are provided in Table 1 and Table S1 of Supplementary Materials.   Table 1 and Table S1 of Supplementary Materials.
The sampling locations were 36 permanent aquatic environments: 9 stagnant water bodies (i.e., lakes including springs, reservoirs and wetlands) and 27 streams ( Figure 1; Table 2). Their altitude ranged from 840 masl to 1600 masl (average slopes of 4-20%). These 36 locations were visited once over a couple of weeks in summer (August 2017). At each of these locations, samples were randomly collected on the different substrates encountered. BDs were sampled by collecting 4 colonized cobbles whose surfaces were around 20 cm 2 ; algal materials were scraped with a brush in a composite sample and preserved in 20 mL of distilled water with 4% formaldehyde. When the substrate was different from that of the paving stones (e.g., submerged macrophyte leaves), a similar colonized surface was sought and BDs were obtained and preserved in the same way. In total, we obtained 132 samples or BD assemblages (Table S1 in Supplementary Materials). The average spacing [33] between the sites in our study was 1.7 km. This discrimination allowed us to carry out analyses both on the total number of assemblages and on sets of them in order to compare the two basins (Júcar vs. Tajo), the different habitats (stagnant water vs. streams) and the different substrates (mineral vs. vegetal). For more details, see Supplementary Materials where more specifications of the studied sites (e.g., sampling location, coordinates, type of habitat and substrate) are presented.

Data Collection Matrices
Independent variables were compiled in four matrices. Sexagesimal latitudes and longitudes were transformed into decimal data resulting in the "geographical matrix" which enabled us to calculate a matrix of Euclidean distances for further statistical analyses. The SigPac Spanish resource (sigpac.mapama.gob.es/fega/visor/, accessed on 21 October 2022) was used to derive data on the altitude, latitude, longitude, area (ha) and mean slope (%) of the sub-catchment upstream of the sampling site. Moreover, because coarser scales of spatial analyses are usually related to catchment features in stream ecology [58], the surface area of main land use types (rangeland and forests) was estimated and is either reported in absolute (ha) or relative (%) terms. Since the territory is highly depopulated (see Supplementary Materials), croplands are now kept to a minimum and they only barely surround some towns. Therefore, cropland areas were included within rangelands. These catchment features comprised the "catchment matrix".
In situ measurements and chemical data obtained in the laboratory encompassed the raw "physico-chemical matrix". Channel cross-section was estimated by planimetry at each sampling station, and the discharge was measured three times at a central point using a FLOW Global Water FP 101 probe. Both measurements enabled us to estimate stream water velocity. Water temperature, dissolved oxygen, pH and conductivity were measured with ODO Yellow Springs and CRISON portable equipment. Chemical analyses using [59] were undertaken shortly after water collection. Nitrogen (nitrate, nitrite and ammonia) and dissolved phosphorus compounds were measured with a Seal analyzer-3, whereas organic carbon (total, TOC hereafter, and dissolved, DOC hereafter) and total nitrogen were measured with TOC-V CSH Shimadzu equipment. Raw samples were digested with strong acids to mineralize all phosphorus forms to render total phosphorus as SRP, which was measured as above. Few data were recorded for silicon, ranging 1.05-2.08 mg Si/L, thus suggesting non-limited diatom growth [60].
The "biological matrix" was constructed with local biological variables which were likely to impinge on BD: presence of riparian arboreal vegetation, presence of herbivory (i.e., the occurrence of fish with a benthic feeding mode), plants as substrate (the details about different specific materials of substrates are provided in Table S1) and evidence of anthropogenic degradation.
BD data were used as dependent matrices. We used a qualitative presence matrix (i.e., presence/absence) to estimate BD richness and frequency of occurrence (i.e., percentage of samples where a taxon was found) and a quantitative matrix with the relative abundance of taxa in each sample. Methods for sample treatment, taxonomical classification and counts were standard and they are explained in the Supplementary Materials file.
Based on these data, statistical analyses could be carried out for the whole or for each catchment (Júcar or Tajo), for each habitat (stagnant water, stream) and for each substrate (mineral, plant).

Population-Community Structure and Its Control
To establish abiotic environmental controls, and suggest either exclusion or coexistence among diatoms, stepwise linear regression models (F-to-enter = 1) were attempted on relative abundance of main taxa and abiotic variables (catchment and physico-chemical matrices) plus biotic variables (the same relative abundance of the main taxa). We defined the main taxa for each sample set (overall, each catchment, each habitat and each substrate) as those with a relative abundance equal to, or greater than, 50% in a given sample, and that were also found in at least 40% of the samples in any of the sets. Before establishing coexistence or exclusion on the basis of a positive or negative relationship between two species, we checked that they did not have a similar relationship with environmental variables. This allowed us to rule out that the latter relationship was the main reason for their exclusion or coexistence. Stepwise linear regression was performed using the STATISTICA 7.0 software package (StatSoft Inc., Tulsa, OK, USA).
We calculated the alpha diversity (richness, Shannon index and effective number) of each local assemblage (considering all taxa). We expressed the Shannon index based on natural logarithms (nats). We transformed it into the effective number [61] which enabled us to perform ANOVAs comparing the average alpha diversity within each set of samples. We also analyzed the abiotic environmental control on alpha diversity with the same approach as above.

Metacommunity Structure, Its Control and Spatial Scales
Based on the composition of the assemblages, we calculated several measurements of βD because their information can be complementary, and it also allows us to compare with results from different authors who used different indices. Based on the presence/absence matrix, Jaccard similarity and Harrison indices of βD [62] were calculated using the PAST package [63]. Multiple-site dissimilarity measurements [31,64], i.e., the spatial turnover (species replacement between sites) and nestedness components of Jaccard dissimilarity, were estimated with the Betadiver R-program (on Jaccard's grounds) in the vegan 2.5-6 package [65]. A nestedness metric based on overlap and decreasing fill (NODF hereafter) was also estimated [64,66]. We also measured βD by taking into account the relative abundance of the taxa [67]: the local (LCBD) and taxa (SCBD) contributions to the quantitative beta diversity. Values were calculated with the beta.div function of the adespatial R-package [68]. Bray-Curtis based ANOSIM was applied to relative abundances of taxa to detect differences in averaged similarity within a set of samples (catchments, habitats and substrates), as mediated by altitude. In order to investigate which taxa contribute the most to the variability in similarity between stagnant and stream waters, we used SIMPER analysis [69]. These analyses were performed using the PAST package [63].
To check the spatial and environmental contributions to BD metacommunity variability, partial RDA (redundancy analysis) for variance partitioning was undertaken following the procedures of [32]. The independent matrices were the spatial (geographical matrix) and environmental ones as reported below (i.e., catchment, physico-chemical and biological matrices). Correlation between environmental variables was calculated to discard collinearity inflation (>80% explained variability to each other). Variance inflation factors were estimated using the vif algorithm from the R vegan package [65]. Finally, the independent components were space and the three environmental matrices separately with their selected variables: physico-chemical (i.e., conductivity and total phosphorus), catchment features (i.e., altitude and the relative area covered by forests) and biological (i.e., substrate type, fish herbivory and riparian canopy). These new matrices were then converted into distance matrices (Euclidean distances) using the varpart program from the vegan package of R [65].
Only those taxa with a relative abundance of more than 3% in more than one sample were included in presence and relative abundance matrices [5]. After being transformed by the Hellinger approach, their Bray-Curtis similarity matrices were calculated and used as dependent matrices in RDA [32]. RDA was separately undertaken for presence and relative abundance as dependent factor.
A complementary analysis to the preceding ones is a spatial correlogram, and this was performed using the adespatial package of R [68]. Spatial autocorrelation was assessed by transforming geographical coordinates into a Euclidean distance matrix. Moran's I and Mantel correlograms were calculated to detect autocorrelation in BD presence and relative abundance matrices. We used the Mantel correlogram function of the vegan R-package [65], and the number of computed distance classes was chosen to follow Sturge's rule [32]. In addition, we also performed the same analysis for an environmental matrix including variables without collinearity (see above). Second-order spatial stationary conditions must, at least, be applied to meet the condition that the spatial variability in data is adequately described by the same single spatial correlation function throughout the study area [32]. This was undertaken by using the resid function of R.

Benthic Diatom Assemblages and Their Controlling Factors
We recorded 156 benthic diatom taxa in the studied area (Table S2). The community richness ranged between 3 and 29 taxa. The average richness was higher in the Tajo catchment compared to the Júcar one, and also in stagnant water compared to streams, but no significant difference was found between epiphytic assemblages and epilithic ones ( Table 3). The range of Shannon index values of all communities was 0.11-2.77 nats, and the range of effective number was 1.12-15.89. The averaged effective numbers between sets of samples did not reveal significant differences (Table 3). Therefore, diversity of taxa was not enhanced by a certain catchment, habitat or substrate. Taking into account all the assemblages, there were no relationships (p > 0.05) between α diversity metrics and abiotic environmental factors. Table 3. Alpha diversity indices, maximal relative abundances and percentage of occurrence (in brackets) of main taxa and beta diversity features of benthic diatom taxa in the water sites within the study area (Serranía de Cuenca) during summer. Results obtained for total samples and on subsets of samples based on catchment, habitat or substrate types. Only species attaining relative abundances greater than 50% in a single community, and occurrences higher than 40% for a given group of samples, are listed as main taxa. Abbreviations: SD: standard deviation; probability (p) of Mann-Whitney test; NODF: nestedness metric based on overlap and decreasing fill. * Not comparable values due to the lower number of sites. Eight taxa made up the most remarkable ones (i.e., dominants and/or cosmopolite; Table 3). Achnanthidium minutissimum, Cocconeis placentula and Gomphonema angustatum were the dominant taxa, with high occurrence in all sets of samples. However, other species seem more exclusive to the catchments, habitats and substrates: Cymbella affinis and Cymbopleura amphicephala mainly occurred in stagnant waters in the Tajo basin, the former mostly appearing epilithically and the latter epiphytically. Cymbella delicatula, Diatoma vulgare and Navicula cryptotenella were mainly observed in streams of the Júcar catchment and on plant substrate.

Overall
The linear models of taxa vs. environmental factors (including all BD taxa plus all variables of catchment, physico-chemical and biological matrices) showed that the variability in BD populations was usually better explained by the relationships with other BDs (i.e., exclusion, coexistence) rather than by any environmental factor (Table 4). For example, in stagnant waters, Cymbella minuta was inversely related to Cymbella affinis and Achnanthidium minutissimum, and it was positively related with Epithemia goeppertiana. These three relationships explained 80% of C. minuta relative abundance. In streams, either A. minutissimum or C. placentula happened to be inversely related with a wide array of other taxa, thus showing some exclusion evidence, albeit the explained variability was around 40% (Table 4). Negative relationships mainly occur between species of the same genus and most positive relationships between species of different genera (Table 4). Only a few populations appear to be somewhat related to the environment, and specifically Cymbella cesatii was strongly controlled by oxygen saturation (>50% of its variability; Table 4). Table 4. Factors (environmental and diatom populations) related to the main taxa of benthic diatom taxa living in the water sites within the study area (Serranía de Cuenca) during summer, as estimated through stepwise multiple regression (p < 0.05). In each section, relationships with abiotic factors are mentioned first, followed by relationships between diatoms.

Metacommunity Structure: Beta Diversity
The large catchments (Júcar and Tajo) shared (Jaccard similarity index) 54% of taxa. Roughly, the same percentage of common taxa was shared between stagnant water bodies and streams. Similarity between communities from different substrate types was analyzed in detail, and the results showed that the filamentous algae substrate shared 47% of taxa with mineral substrate and 48% with macrophytes, and that mineral substrate and macrophytes shared 57% of taxa. When split altitudinally by 200 m ranges, Jaccard index values were also between 47% and 53%.
The Harrison βD was similar in both catchments, whereas that of stagnant waters and mineral substrates were greater than that of streams and epiphytic substrates, respectively ( Table 3). The range of the NODF metric was 22.7-29.7, with stagnant water bodies showing the highest value ( Table 3). The spatial turnover range was 0.73-0.96, and dissimilarity due to nestedness varied between 0.02 and 0.11 (Table 3). The SCBD and LCBD ranges were 0.000-00.123 and 0.005-00.013, respectively. Only the species contribution (SCBD) to βD was significant.
Considering the relative abundances of the taxa, the averaged similarity (ANOSIM, Bray-Curtis index) was slightly higher in the Tajo compared to the Júcar catchment (0.26 ± 0.20 vs. 0.21 ± 0.22; R 2 = 0.13, p = 0.0001) and in stagnant waters compared to streams (0.28 ± 0.25 vs. 0.22 ± 0.20; R 2 = 0.13, p = 0.02). Differences between substrates were not significant (p = 0.09). There was an altitudinal effect, as at intermediate altitudes (1000-1400 masl) the averaged similarities were lower than those averaged at higher and lower altitude sites (Figure 2). To detect a possible effect of the catchment on this altitudinal pattern, we compared the relative abundance of the eight distributions (two catchments x four altitudinal ranges; Figure 2) which resulted in significant differences (Kruskal-Wallis test: H-Chi 2 = 107, p < 0.0001). The Mann-Whitney post hoc pairwise comparisons showed that the mean similarity at intermediate altitudes, regardless of catchment, was lower than at higher altitudes ( Figure 2). Seven taxa contributed up to 60% of the dissimilarity between stagnant water and streams in each catchment (Table S3), four of which were common to both catchments: Achnanthidium minutissimum, Cocconeis placentula, Cymbopleura amphicephala and Gomphonema angustatum.
The Harrison βD was similar in both catchments, whereas that of stagnant waters and mineral substrates were greater than that of streams and epiphytic substrates, respectively ( Table 3). The range of the NODF metric was 22.7-29.7, with stagnant water bodies showing the highest value ( Table 3). The spatial turnover range was 0.73-0.96, and dissimilarity due to nestedness varied between 0.02 and 0.11 (Table 3). The SCBD and LCBD ranges were 0.000-00.123 and 0.005-00.013, respectively. Only the species contribution (SCBD) to βD was significant.
Considering the relative abundances of the taxa, the averaged similarity (ANOSIM, Bray-Curtis index) was slightly higher in the Tajo compared to the Júcar catchment (0.26 ± 0.20 vs. 0.21 ± 0.22; R 2 = 0.13, p = 0.0001) and in stagnant waters compared to streams (0.28 ± 0.25 vs. 0.22 ± 0.20; R 2 = 0.13, p = 0.02). Differences between substrates were not significant (p = 0.09). There was an altitudinal effect, as at intermediate altitudes (1000-1400 masl) the averaged similarities were lower than those averaged at higher and lower altitude sites (Figure 2). To detect a possible effect of the catchment on this altitudinal pattern, we compared the relative abundance of the eight distributions (two catchments x four altitudinal ranges; Figure 2) which resulted in significant differences (Kruskal-Wallis test: H-Chi 2 = 107, p < 0.0001). The Mann-Whitney post hoc pairwise comparisons showed that the mean similarity at intermediate altitudes, regardless of catchment, was lower than at higher altitudes ( Figure 2). Seven taxa contributed up to 60% of the dissimilarity between stagnant water and streams in each catchment (Table S3), four of which were common to both catchments: Achnanthidium minutissimum, Cocconeis placentula, Cymbopleura amphicephala and Gomphonema angustatum.

The Metacommunity: Variance Partition and Spatial Patterns
The RDA showed that both spatial and environmental factors explain both the presence and the relative abundance of taxa, and do so with similar efficiency (Table S4). We will, therefore, only detail the results on the relative abundance matrix. Unexplained variance was around 50% (Table 5 and Table S4). The pure environmental component explained at least twice the variance of relative abundance compared to the pure space component (Table 5); the rest was due to the overlap of space with environmental contributions (Figure 3a). When using the three environmental matrices (physico-chemical, catchment and biological matrices), it was the biological matrix (i.e., substrate type, fish herbivory and riparian canopy) that explained some variability independently of space (Figure 3a; Table 5), and pure space did not explain any variance. The relative abundance data, grouped by catchments, habitat or substrate (Table 5), showed less spatial intervention and slightly more variance explained in streams compared to stagnant waters and in the Júcar catchment (with a higher proportion of streams) than in the Tajo, as well as in the epilithic samples for which the highest values of variance explained by the pure biological component were observed. Table 5. Spatial and environmental relative contribution (%) to total explained variation of relative abundance of benthic diatom taxa. RDA results correspond to overall samples and different sets of samples and performed on spatial and environmental factors or on spatial and three environmental matrices (physico-chemical, catchment and biological). Pure means without interaction with other matrices. Only statistically significant results (p < 0.05) are shown.  Figure 3. (a) Spatial and environmental (catchment, physico-chemical and biological matrices) contributions to BD metacommunity variability, calculated on relative abundance of taxa. Unexplained variance was 51%; the distribution of the variance explained is plotted. (b) Moran's I spatial autocorrelation coefficients of the environmental and the relative abundance of BD taxa. Eight distance geographical classes (from 0 to 46 km, with each class being 6 km) were analyzed; only statistically significant coefficients at p < 0.05 are plotted. Figure 3. (a) Spatial and environmental (catchment, physico-chemical and biological matrices) contributions to BD metacommunity variability, calculated on relative abundance of taxa. Unexplained variance was 51%; the distribution of the variance explained is plotted. (b) Moran's I spatial autocorrelation coefficients of the environmental and the relative abundance of BD taxa. Eight distance geographical classes (from 0 to 46 km, with each class being 6 km) were analyzed; only statistically significant coefficients at p < 0.05 are plotted.

Benthic Diatom Populations-Communities and Their Control Factors
The geographic variation of the BD community (Mantel's spatial autocorrelation with Moran's I coefficient across distance classes) was considerably different whether calculated with the species occurrence matrix or the relative species abundance matrix (Table S5). BD presence only showed autocorrelation (p = 0.028) at the shortest distance classes, and BD relative abundances followed a decreasing curve (Table S5, Figure 3b). When we estimated the autocorrelation pattern of environmental factors used above, a U-shaped pattern appeared (Table S5, Figure 3b).

Benthic Diatom Populations-Communities and Their Control Factors
The paucity of BD studies has prevented the species-area relationship from being confirmed so far with a meta-analysis [47]. In fact, the disparity of these values can be seen in Table 1. For example, the richness reported here (with a study area of 3300 km 2 ) is similar to that found in the large area of the Manyame catchment in Zimbabwe (40,000 km 2 ; [14]) or in a smaller region in China (950 km 2 ; [15]). To this hypothetical effect of the area, we should add the effect of the type of habitat that is mainly included in such areas. For example, [53] and ourselves, studying more than twice as many samples from streams as from stagnant water, found 25% richer communities in the latter habitats than in the former. In addition, the Júcar catchment, whose rivers:lakes ratio of samples is three times greater than that of the Tajo catchment, has the lowest richness in its communities. These results would agree with the expected lower richness in highly disturbed locations [37,70].
Furthermore, we wonder whether the type of substrate which is usually sampled should also be taken into account. So far, there is much more data on rock communities than on plants or other mineral substrates (Table 1). Perhaps subtler differences will be found when more data on substrates become available. For the time being, the richness of communities on plants and on minerals according to our results was similar, and this is also inferred from other works (e.g., [71][72][73]). In addition, we have proved that ecological diversity, measured as the effective number of species, was independent of the habitat and substrate, as other studies have also reported for the Shannon index [53,72,73]. Therefore, the fact that the selection of the sampled substrate is certainly biased towards cobbles (Table 1) does not seem to be relevant in establishing the ecological diversity of BD. At the moment, it does not seem to be different from the diversity of assemblages from, for example, plants. This is noteworthy because, otherwise, the effort to obtain the ecological diversity of assemblages from various substrates of an ecosystem would have to be extended, increasing the costs of, for example, water quality assessment. Furthermore, in the intermediate spatial extent studied here, alpha diversity does not appear to reflect changes in any environmental factors. Our first hypothesis, "α diversity will be independent of environmental factors and favored in undisturbed systems", has been confirmed. In addition, we discovered as a corollary that the variety of habitats, and the intensity of disturbances that are amplified with the area studied, will make it difficult to find a BD species-area relationship.
Our results suggest that interactions between BD populations at the local scale should be taken into account as mechanisms involved in BD assemblage construction. These biotic relationships explain the variability in taxa distribution better than the abiotic environment. One of the main results regarding likely mechanisms implying the exclusion effect (or overdispersed distributions in the sense of [37]) was that of the opposite dominances of Achnanthidium minutissimum and Cocconeis placentula in streams, something which had already been mentioned in the 1970s [40]. More recently, [42] attributed the contrasting distribution of these two populations to a greater dispersal ability of the small A. minutissimum, a mass effect under environmentally sub-optimal conditions, versus C. placentula, which is a larger diatom and is more prone to a clumped distribution. These two dispersive strategies based on their traits allow these species to also be ubiquitous, appearing in different continents and dominant in small areas, such as the one analyzed here, or in other much larger ones ( [4]; Table 1). Another mechanism would explain the contrasting distribution of Cymbella affinis and Cymbella minuta in stagnant waters (Table 3), namely, competition between very similar and fast-growing species [74], such as these ones. It has been impossible to find studies on the biology of these Cymbella species that can be related to their overdispersion distributions due to competition. Despite the fact that competition between benthic algae, as a mechanism that would explain their distribution, was suggested at the end of the 1990s [41], interactive mechanisms between BD have been scarcely addressed [43,75]. We have found positive relationships between some taxa such as Cymbella minuta and Epithemia goeppertiana (Table 3); this is possible because these two taxa have very different strategies to take advantage of the low phosphorus concentration [50,76]. In addition, the fact that most of the positive relationships occur between species of different genera, and the negative ones mainly between species of the same genus, would reinforce the idea that mechanisms avoid sharing the niche work. We consider our second hypothesis, "relationships between BD populations will become very important for the assembly process, with exclusion (overdispersion) being the most relevant", to have been tested. Moreover, we can express another corollary: there is a small set of dominant and ubiquitous species, but they belong to different assemblages as is evident from their preferences for rivers or lakes and their overdispersion distribution in communities. This assemblage design, the coexistence of taxa recorded in the different BD communities, is reminiscent of the species associations recognized for phytoplankton which were related to the environment [77,78]. This would be an interesting avenue to pursue with BD.

Metacommunity Structure: βD Patterns
The shared BD taxa among different groups of samples was around 50%. This ubiquity of BD was independent to catchments, habitats or substrates. This component of homogeneity in communities results in a low Harrison index and weak nestedness, both indices being lower in streams than in stagnant waters. Therefore, we suggest that at the intermediate scale (two contiguous basins) studied here, it is possible to observe how the water flow homogenizes the environment that BDs inhabit. This result, along with those of [24,79], point to how homogeneous pristine rivers seem to be for diatom distribution when analyzed at these sparsely studied intermediate scales. In addition, a higher homogeneity in BD communities in streams compared to stagnant waters has been observed on other occasions (e.g., [53]), and the cause may be that higher physical connectivity in streams favors dispersion of microorganisms and environment homogenization (e.g., [9,80]).
On the other hand, a high spatial turnover and a negligible value of the nestedness component of dissimilarity means that the βD pattern in our intermediate studied region was almost exclusively caused by species replacement [64]. The spatial turnover component was greater in our study than the value recorded comparing numerous French streams at a larger scale [23]. In addition, the use of relative abundance allows us to highlight that the contribution of species to βD (SBDC) is considerably higher than the local ecological uniqueness contribution (LCBD), the latter being related to the unusual habitat conditions that select only a few taxa [81]. While the range of SCBD at the intermediate spatial scale used here was greater than that of other studies of diatoms over a large regional scale, our range of LCBD was within their reported ranges [19,26,29]. The SCBD reflects the fact, already described here, that there is a small group of dominant species (only seven species are responsible for the dissimilarity) which are also the most ubiquitous and are exclusive of each other (e.g., A. minutissimum and C. placentula), and all other species are rare. This seems to be the BD metacommunity pattern, both at intermediate and large spatial scales such as the one studied by [19]. Therefore, the mechanisms underlying this pattern, supported by the above analysis of population controlling factors and compositions of communities, seems to be biotic interactions between BDs in the assembly process instead of a selective environment factor acting as a filter for main populations. The assembly process as a mechanism of main species sorting is reflected in the high turnover and rare species following a mass effect or a stochastic distribution. As the replacement pattern predominates between BD communities, the maintenance of their biodiversity involves the conservation of communities of different localities with a criterion of floristic complementarity [26,31].
When looking at βD at different altitudes, we found other possible evidence that higher environmental stability implies greater beta diversity as a result of the assembly process: the highest βD occurs among the communities at intermediate altitudes. Despite the reduced altitudinal range (840-1600 masl) studied here, the climatic gradient could be considered a driver of altitudinal variation in βD [47] because temperature was one of the few local environmental factors explaining variance in some BD populations. At headwaters, or higher altitude, climate acts as a filter, while downstream, or at lower altitudes, high pollution or anthropic degradation would be another filter (nitrogen being the other somewhat significant factor for populations).
Therefore, we consider our third hypothesis (i.e., βD will be higher if environmental stability favors alternative assemblages), linking stability to the possibility that assemblage processes lead to different communities by combining a set of species, to be proven. Even in a small-intermediate area with little altitudinal variation, we detected lower βD in lotic habitats compared to lentic habitats and at more disturbed altitudes. Further studies on BD are needed to understand the relationship between all the different βD metrics and their underlying explanatory mechanisms and processes [9,43,[82][83][84]. For the time being, to unravel the βD components [67] of BD metacommunities the use of pre-existing databases could be of interest (Table 1).

The Metacommunity: Variance Partition and Spatial Patterns
Presence and relative abundance were explained by both space and environment, as reported by other studies (Table 1). Relatively high unexplained variability often occurs in many metacommunity studies as a result of overlooked factors, stochasticity, interactions with unconsidered taxa, etc. [26,85]. Presence and relative abundance both showed similar good results (Table S4) and are therefore recommended for the monitoring of BD metacommunities [26,47]. Pure environment, without space overlap, explained double the variance compared to pure space (Table 5; Figure 3a). The importance of environmental factors for diatom metacommunities has already been emphasized by [10]. Ref. [12] suggested that it is in small areas, such as those used in this study, where more variability purely explained by the environment is experienced, because geographical features in larger extents may mask it. Related to this, pure space did not explain the variance, or explained very little variance in the Serranía of Cuenca area. In a slightly larger area, pure space was found to be significant ( [5]; Table 1), but the author suggested that this spatial variance was due to a spatial gradient of land use attributes. In our study, landscape or catchment variables, such as land use, were already included in environmental matrices, and therefore our result concerning structuring due strictly to space clearly shows the dispersive capabilities of BD populations. To our knowledge (Table 1), this is the smallest extent to which pure space explained some, albeit little, variance in BD. We verified that it is in streams where metacommunity structure was better explained, especially by the environment, as was also observed by [10]. When treating space together with physico-chemical, catchment and biological matrices separately, the biological variables (i.e., riparian canopy, substrate types and fish herbivory) were the ones that purely explained the variance. This fact highlights that biological relationships at the local scale are important in the structuring of the metacommunity [10]. All these results encourage the inclusion of more biotic variables, such as the presence of bacteria, other benthic algae or macroinvertebrates strongly linked to BD growth [43,86] or the three-dimensional mode of their community [87]. We should test whether the addition of new local biotic variables, relevant to the growth of BD, help us to understand the drivers in the structuring of their metacommunity.
The weak effect of space versus environment suggests that the BD metacommunities studied here arise largely from species sorting, and from species sorting plus mass effect processes [30,34,85], as inferred from the above βD analysis and as expected in small areas [47], confirming our fourth hypothesis. In conclusion, we found that local interspecific relationships, or the process of community assembly, were the mechanisms explaining βD (i.e., we highlight the relevance of interspecific biotic relationships) and that BD relationships with local biological conditions shape the metacommunity structure and explain its construction pattern.
In addition, the above conclusions on the relevance of local factors structuring communities and the metacommunity was reinforced when disentangling the patch structure over different spatial scales. An initial positive value of autocorrelation at the smallest scale (less than six kilometers) was followed by significant negative values at larger scales, tens of kilometers, as detected by [6] in a fragmented French river. The correlogram of the environment also showed similar positive values at the finest spatial scale, confirming that the environment structures the spatial patterns of the BD at this scale (zone of influence, sensu [88]).
Finally, we would like to draw attention to the fact that both the metacommunity archetypes (species sorting, species sorting plus mass effect) and the scale of interaction with the environment have been described not only for small but also for large extents and coarser scales [4,16] addressing similar goals about the spatial ecology of BD (Table 1). This may be a serious inconsistency that has long been detected in both geography and ecology [89,90] and may prevent cross-comparison of studies [33]. In short, we need to understand the environmental factors acting at the microscopic scale as determinants of organismal ecology (e.g., inside the biofilm; [43]), the local factors of the environment where BD assembly occurs and the community develops [57] and at a broader scale the environmental patterns defining macroecological rules [91]. Based on the confirmation of our fourth hypothesis, we encourage the analysis of the local communities' drivers to understand the metacommunity structure.

Conclusions
As in an experimental design in which almost all conditions are set to analyze the variability in a variable in depth, we wanted to show, with the use of an intermediate spatial extent and the inclusion of many biotic variables, that biotic relationships occur that explain the presence of populations and the organisation of communities that ultimately structure the metacommunity. As a result of the implementation of conservation directives and concerns about global change, there are a large number of quantitative BD datasets available worldwide that could be a very good basis for (following our recommendation) achieving a better understanding of BD ecology. Some conclusions of our work, due to their relevance in applied ecology, should be tested in more environments, e.g., the indifference of assembled BD to different substrates and the poor response of diversity indices to environmental changes. Finally, we call for future work on the analysis of the relationships that structure communities (populations with their biotic and abiotic environment) at local scales as a basis for understanding the mechanisms underlying environmentally explained variance in metacommunity structure.