Landscapes of facilitation: how self-organized patchiness of aquatic macrophytes promotes diversity in streams

Spatial heterogeneity plays a crucial role in the coexistence of species. Despite recognition of the importance of self-organization in creating environmental heterogeneity in otherwise uniform landscapes, the effects of such self-organized pattern formation in promoting coexistence through facilitation are still unknown. In this study, we investigated the effects of pattern formation on species interactions and community spatial structure in ecosystems with limited underlying environmental heterogeneity, using self-organized patchiness of the aquatic macrophyte Callitriche platycarpa in streams as a model system. Our theoretical model predicted that pattern formation in aquatic vegetation – due to feedback interactions between plant growth, water flow and sedimentation processes – could promote species coexistence, by creating heterogeneous flow conditions inside and around the plant patches. The spatial plant patterns predicted by our model agreed with field observations at the reach scale in naturally vegetated rivers, where we found a significant spatial aggregation of two macrophyte species around C. platycarpa. Field transplantation experiments showed that C. platycarpa had a positive effect on the growth of both beneficiary species, and the intensity of this facilitative effect was correlated with the heterogeneous hydrodynamic conditions created within and around C. platycarpa patches. Our results emphasize the importance of self-organized patchiness in promoting species coexistence by creating a landscape of facilitation, where new niches and facilitative effects arise in different locations. Understanding the interplay between competition and facilitation is therefore essential for successful management of biodiversity in many ecosystems.


INTRODUCTION
The challenge of understanding species diversity and coexistence is fundamental in community ecology. According to the competitive exclusion principle, two species competing for the same resource cannot coexist if other ecological factors are constant (Gause 1932). However, many natural communities defy the theoretical predictions of low species coexistence, as often a high number of species can be found living on few resources (e.g., "paradox of the plankton", Hutchinson (1961)). To explain this discrepancy, many of the suggested mechanisms rely on the importance of spatial or temporal heterogeneity (Levin 1970, Koch 1974, Armstrong and McGehee 1976, Holt 1984, Tilman 1994, Amarasekare 2003. Extensive evidence exists that structurally complex physical habitats favour increased species diversity, by providing niches and diverse ways of exploiting environmental resources (MacArthur and MacArthur 1961). Yet, many ecosystems with limited abiotic heterogeneity also host a high number of species. Thus, despite the importance of heterogeneity in space and time for species coexistence, we still lack understanding of how species can coexist in environments where underlying abiotic heterogeneity is low.
In recent decades, there has been increasing evidence that strong interactions between organisms and their environment can create environmental heterogeneity, even under uniform, homogeneous conditions, through the process called spatial self-organization (Sol e and Bascompte 2006, Rietkerk and Van de Koppel 2008). Self-organization processes can generate spatial patterns in ecosystems, through the interaction between local positive and large-scale negative feedbacks (Rietkerk and Van de Koppel 2008). Examples range from vegetation patches alternating with bare soil areas in arid ecosystems (Rietkerk et al. 2002), tree patterns in Siberian peatlands (Eppinga et al. 2008) to diatoms in homogeneous tidal flats (Weerman et al. 2010). Self-organized patterns can cause strong variability in abiotic conditions in their surroundings. By modifying the abiotic environment, self-organizing species can promote favourable conditions leading to a positive feedback on their own growth (Wilson and Agnew 1992, Rietkerk and Van de Koppel 2008, K efi et al. 2016.
Several studies have also focused on the importance of positive interactions that benefit individuals of different species, i.e., interspecific facilitation (Bertness and Callaway 1994, Pugnaire et al. 1996, Callaway and Walker 1997, Brooker et al. 2008. For instance, facilitator species can reduce environmental stress, increasing the realized niche of other species and allowing them to occupy environments that they would normally not inhabit (Bruno et al. 2003, Callaway 2007. Facilitation is in essence based on the same mechanism as self-organization, involving a positive interaction that improves environmental conditions and enhances growth or survival. However, facilitative interactions between two species are mostly considered at a relative local scale, within a tussock or patch of the facilitator species, for instance through "nurse plant effects" in relation to herbivory or drought (Callaway 1995, Padilla andPugnaire 2006). Instead, studies of self-organization typically focus on a single species at a landscape setting, analysing both scaledependent effects of local facilitation and large-scale competition (Rietkerk and Van de Koppel 2008, van Wesenbeeck et al. 2008, Schoelynck et al. 2012. Therefore, as the link between self-organization and interspecific facilitation remains unclear, we pose the question whether self-organized pattern formation can create a "landscape of facilitation". In lotic aquatic ecosystems, self-organized patchiness has been found to occur in submerged aquatic vegetation due to scale-dependent feedbacks between plant growth, water flow and sedimentation processes (Schoelynck et al. 2012(Schoelynck et al. , 2013. Submerged macrophytes often grow as well-defined, streamlined stands composed of either a single species, or a mixture of species. Macrophytes act as ecosystem engineers (Jones et al. 1994), slowing down the water flow within the patches and promoting sediment deposition (Sand-Jensen and Mebus 1996, Sand-Jensen 1998, Wharton et al. 2006, which creates a local positive feedback on their own growth and survival. At the same time, flow velocities increase around the patches, creating a large-scale negative feedback on plant growth due to the increased mechanical stress (Puijalon et al. 2011, Schoelynck et al. 2012. In lowland rivers, aquatic macrophytes with different morphologies increase habitat heterogeneity beyond that promoted by hydrodynamic and geomorphological processes alone (Kemp et al. 2000, Gurnell et al. 2006. Despite being suggested by previous observational studies (Jones 1955, Haslam 1978, the consequences of such plant-driven heterogeneity for interspecific interactions have not yet been explored. We investigated whether self-organized pattern formation in aquatic vegetation promotes the coexistence of different macrophyte species in lotic communities, by generating heterogeneous hydrodynamic conditions and hence creating a "landscape of facilitation". First, to demonstrate self-organized pattern formation by the aquatic macrophyte Callitriche platycarpa K€ utz (various-leaved water starwort), we constructed a spatially explicit mathematical model based on the interaction between plant growth and hydrodynamics. Secondly, we investigated whether such self-organized spatial heterogeneity could promote species coexistence, by modelling the interaction between the pattern-forming species (i.e., facilitator) and two species (i.e., beneficiaries) with different resistance to hydrodynamic stress. Thirdly, to show self-organization and spatial association among species in the field, we compared the model-predicted spatial distribution patterns against field observations on the spatial distribution of two hypothesized beneficiary species (lesser water parsnip, Berula erecta (Huds.) Coville and opposite-leaved pondweed, Groenlandia densa (L.) Fourr.) around Callitriche. Finally, to show that such spatial association provides facilitative interactions, we carried out field transplantations of the two beneficiary species in different locations around patches of the facilitator Callitriche as well as on bare sediment, and we investigated if their growth rate, reproduction, and survival correlated with changes in hydrodynamic conditions created by Callitriche patches. Our results suggest that species coexistence in streams is promoted by a biophysical feedback process that creates a landscape of facilitation where multiple new niches emerge for species adapted to a wide range of conditions.

MATERIALS AND METHODS
A model of pattern formation for submerged aquatic macrophytes Model description.-To study the emergence of self-organized patterns in aquatic macrophytes and the potential consequences for species coexistence, we constructed a spatiallyexplicit mathematical model based on the feedback between vegetation and water flow. The model consists of a set of partial differential equations, where two equations describe the dynamics of plant biomass for the facilitator species f (P f ) and for its beneficiary species b (P b ), and where water velocity in the streamwise and spanwise directions (u and v), and water depth (h) are described using the shallow water equations (Vreugdenhil 1989). The rate of change of plant biomass per species in each grid cell can be expressed as: where i = f and j = b for the equation of the facilitator (pattern-forming) species, and vice versa for a beneficiary (nonpattern forming) species. Here plant growth is described using the logistic growth equation, where r i is the intrinsic growth rate of the plants and k i is the plant carrying capacity.
Competitive interactions between P f and P b are accounted for using the competitive Lotka-Volterra equations, with the term a ij representing the effect P j has on P i . Plant growth rate r i is reduced when sediment accumulation within the plants increases towards its maximum value A; this represents a negative feedback on plant growth due to sediment accumulation and organic matter content becoming high enough to be toxic for the plants ( Barko and Smart (1983); Sofia Licci, personal communication). S is the sediment level (m). Plant mortality m is assumed to decrease with increasing plant density because of a reduction of flow stress in dense vegetation. This is represented by the term F/(P i + F), where F is an intraspecific facilitation term. Plant mortality caused by water flow stress is modelled as the product of the mortality constant m Wi and net water speed u j j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu 2 þ v 2 Þ p (m/s) due to plant breakage or uprooting at higher velocities (where u and v are water velocities in the streamwise and spanwise directions). Field sampling on clonal dispersal traits for the aquatic plant species Berula erecta and Groenlandia densa revealed that plant lateral expansion through vegetative reproduction could be described by a random walk (see Appendix S1: Fig. S1). Therefore, we apply a April 2018 SELF-ORGANIZED LANDSCAPE OF FACILITATION 833 diffusion approximation and use these data to parameterize different diffusion constants D j for the beneficiary species (Holmes et al. 1994). Clonal dispersal traits for the hypothesized facilitator species Callitriche platycarpa could not be estimated based on field sampling, due to the complex morphology of this species. Therefore, the diffusion constant D i for the facilitator species was given an estimate value. Changes in sediment level are described as: where S in is the sediment deposition rate (m/t), E max is the maximal erosion rate of sediment (per t) and K S represents the effects of plants in promoting sediment deposition. The term u j jrS represents the advective flux of sediment over the bottom (i.e., as fluid mud) in any horizontal dimension, and D S represents the horizontal dispersion rate of sediment, mainly due to flow heterogeneity, and to a lesser extent processes such as bioturbation, which is modelled with a diffusion approximation.
Water flow is modelled using depth-averaged shallow water equations in non-conservative form (Vreugdenhil 1989), to determine water depth and its speed in both x and y directions (see Appendix S2 for the complete set of equations and description of the variables). The effects of bed and vegetative roughness on flow velocity are represented by determining hydrodynamic roughness characteristics for each cover type separately using the Ch ezy coefficient, following the approach of Straatsma and Baptist (2008) and Verschoren et al. (2016).
Within the unvegetated cells of the simulated grid, the Ch ezy roughness of the bed (C b ) is calculated using Manning's roughness coefficient through the following relation: where n is Manning's roughness coefficient for an unvegetated gravel bed channel (s/[m 1/3 ]) and h is water depth (m). For each grid cell occupied by submerged vegetation, C d is calculated using of the equation of Baptist et al. (2007) and slightly modified by Verschoren et al. (2016) to account for reconfiguration of flexible submerged macrophytes, to express vegetation resistance as: where C b is the Ch ezy roughness of the bed, g is acceleration due to gravity (9.81 m/s 2 ), D c is a species-dependent drag coefficient, A w is the specific plant surface area ( Table 1 provides an overview of the parameter values used, their interpretations, units and sources. Model analysis: simulation of species coexistence patterns.-To investigate whether spatial pattern formation could promote species coexistence through the creation of spatial heterogeneity in hydrodynamic conditions, we modelled the  Ecology,Vol. 99,No. 4 interaction between the pattern-forming species (P f ; facilitator) and two non-pattern forming species (P b , beneficiary species). In the first model, we considered the interaction between P f and a beneficiary species P b1 characterized by low resistance to hydrodynamic stress (= high mortality constant m Wi ; Table 1). In the second model, we considered the interaction between P f and a second beneficiary species P b2 characterized by higher resistance to hydrodynamic stress (= low mortality constant m Wi ), but lower growth rate and lower dispersal ability. We modelled the pairwise interactions between the facilitator and each beneficiary separately instead of with a full three-species model. This choice was made to focus on the mechanisms and patterns allowing the coexistence of single beneficiary species with the self-organizing species, instead of studying the coexistence patterns of a whole community. Hence, we focused on studying a self-organized landscape with spatial facilitation, rather than exploring all possible modes of coexistence. The models were analyzed by simulating the spatial development of vegetation after random seeding (increasing biomass to 1 in randomly chosen cells) on a spatial grid of 300 9 60 cells, corresponding to a river stretch of 25 m 9 5 m. We investigated vegetation development with two-dimensional numerical simulations using the central difference scheme on the finite difference equations. The simulated area consisted of a straight channel with rectangular cross-sectional shape and initial bed slope of 0.03 m/m. Simulations were started by specifying an initial value of inflowing water speed for the streamwise water flow in the x direction and assuming constant flux. The model was implemented in Matlab (version 2016b, The MathWorks, Inc.). Simulations were run for 500 time steps, in abstract units due to our non-dimensional description of plant growth.
To test the regularity of the predicted spatial patterns, we analyzed the resulting distribution patterns of P f through spatial autocorrelation. To test the spatial dependence between the beneficiary species P b and P f , we used spatial cross-correlation. Both auto-and cross-correlation analyses were performed by calculating Moran's I in the "ncf" package in R (Bjornstad, 2016). To test for self-organization and spatial association among species in the field, we then compared the auto-and cross-correlation functions from the predicted species distribution patterns of coexistence with field observations on the spatial distribution of Groenlandia and Berula around Callitriche (see paragraph 2.2).
To further explore the implications of self-organization for species coexistence, as opposed to homogeneous environments, we compared the spatial model described above to a simplified, homogeneous (non-spatial) version of the model based on Eq. 1: where i = f and j = b for the equation of the facilitator species, and vice versa for a beneficiary species. We used the model to explore the realized niche of each species along the hydrological gradient, under homogeneous (non-spatial) conditions (that is, without self-organization). This simplified version of the model does not account for spatial effects of sedimentation or velocity and intraspecific facilitation.
For each imposed flow velocity U in (|u| in Eq. 5), we explored the conditions under which the model predicted either stable coexistence, unstable coexistence or competitive exclusion between the facilitator and beneficiary species (based on the species isoclines of zero growth), as a result of their stress resistance and competitive abilities. Moreover, to show the hydrodynamic heterogeneity generated by the selforganization process and the species hydrological niches predicted in the spatial model, we investigated the frequency distribution of flow velocities within vegetated and unvegetated cells in the spatial model. The comparison between the two models provided insight and understanding of the mechanisms underlying species coexistence in space.

Field observation of species coexistence patterns through aerial photographs
To test for significant spatial association of species around self-organized patterns in the field, we examined the distribution of two potential beneficiary species (Groenlandia and Berula) around the hypothesized facilitator species (Callitriche). Submerged macrophytes often grow as well-defined stands composed of a single species or a mixture of species (Fig. 1A); the patches tend to merge into a more homogeneous cover where streams have low flow velocities sustained over time (Fig. 1B), while distinct streamlined patches are usually found in streams with sustained periods of moderate to high flow velocities (Fig. 1C). Vegetation distribution was mapped in two reaches of 100 m in length, through lowaltitude aerial photographs. The channels are located along the Rhône River (France), near Serri eres-de-Briord (45.815311°N, 5.427477°E) and Fl evieu (45.766738°N, 5.479622°E) (see Appendix S3: Fig. S1 for the location of the study sites). The first reach was mainly colonized by Callitriche and Groenlandia, with few patches of other macrophyte species, while the second reach was colonized only by Callitriche and Berula. Aerial pictures of the streambed were taken with a digital camera mounted on a pole at about 2 m height that was moved in the upstream direction along the stretch. Aerial pictures were collected at times of day when the sun was at its highest point, and in the few hours before and after it (between 10:00 and 15:00 h), to minimize glare. Pictures were collected with a slight overlap and afterwards mosaicked using image processing software (Adobe Photoshop CC 2015). Patches of different species were identified and delineated as shown in Fig. 1A; afterwards, pixels where the species was absent were given a value of 0 and pixels where the species was present were given the value of its blue channel in the RGB image, since the intensity of this channel was the one most closely related to differences in plant biomass (evaluated by visual inspection). This allowed us to obtain different raster maps of macrophyte distribution, one for each of the species considered in the study (non-target species were not included in the analysis). The resulting macrophyte maps were analyzed through spatial autocorrelation (to test the distribution of the potential facilitator species) and cross-correlation (to test the spatial dependence between the facilitator and each of the potential beneficiary species), by calculating Moran's I. The sample size of our field observations was constrained by the time-intensive nature of our image collection method. Field studies of this nature are often constrained in terms of sample size but can provide valuable insights even without replication (Colegrave and Ruxton 2017). By integrating multiple approaches, the aim of our study was to provide a "proof of principle" for the mechanisms underlying self-organization and species coexistence in space. Furthermore, the field study in a simplified channel provides a valuable starting point for more observations with different aquatic species and in different stream types. Throughout this paper, the term wake is used to indicate a region of reduced velocity directly downstream from a vegetation patch, i.e., where the flow is laterally uniform and slower than the flow around the patch (Zong and Nepf 2012, Liu and Nepf 2016).

Testing for positive interactions through a field transplantation experiment
To test for the presence of positive interactions between the hypothesized facilitator C. platycarpa and the two hypothesized beneficiary species living in its surroundings, we performed a field transplantation experiment in an artificial drainage channel with natural colonization by aquatic vegetation. The channel is located along the Upper Rhône River (France), near Serri eres-de-Briord (45.810657°N, 5.447169°E), it is 4.26 km long, uniform in terms of width and water depth, with relatively straight banks. The average width is 8.0 m and the average depth is 0.8 m, rarely exceeding 1.3 m. The channel has a substrate of fine sand (d 50 = 230.87 lm). Flow velocities are on average 0.25 m/s, with a discharge of 1.48 AE 0.022 m 3 /s as measured on 21 August 2014 (averaged over five transects in the study site). The channel is fed by groundwater supply (see description of the flow conditions in the paragraph after the next one and Appendix S4).
Individuals of the two beneficiary species were collected within the same channel on 11th August 2014 and transplanted in five locations around the facilitator patches. Along the patch central axis, transplants were located 20 cm upstream of the leading edge, in the middle (50% of the patch length) and 20 cm downstream of its rear edge. Next to the patch, transplants were positioned at 20 cm to the left and to the right side of its lateral edges, at 50% of the patch length. As a control, an additional treatment was located on bare sediment areas, as far as possible from the influence of existing patches. Since patch effects can be observed for a distance equal to its length (Sand-Jensen and Mebus 1996, Schoelynck et al. 2012), these transplants were located at a distance of at least twice the length of the nearest patch. Ten transplants per treatment were used for each beneficiary species, with one transplant per position around different C. platycarpa patches of average length (~1.2 m) and in areas outside the influence of other vegetation. Transplants were single plants attached to a stolon without internodes (shoot height of 22.17 AE 1.98 cm for B. erecta, 21.48 AE 1.98 cm for G. densa). All field transplantation experiments cause a disturbance to the system. However, this disruption and thus its impact on the subsequent observations and measurements was kept to a minimum by creating a small hole of approximately 8 cm in depth and 2-3 cm in diameter in the sediment using a metal pole, to accommodate the rooting part of each single plant shoot. The hole was refilled with sediment almost immediately, and small cobbles were placed around it to prevent scouring and washout of the planted shoots. We observed a limited release of sediment at the time of planting and the conditions stabilized within 2 d. Transplant survival was monitored 2 d, 4 d, and at weekly intervals after transplantation to test for facilitative effects on plant survival. All transplanted individuals were harvested at the end of the experiment (49 d after transplantation, on 29th September 2014). The duration and timing of the experiment were designed for a period long enough to enable transplants to grow and reproduce by clonal growth (Puijalon et al. 2008, Schoelynck et al. 2012, and to harvest plants at the end of the growing season, before autumnal decay. No storms took place during the experimental period.

836
LORETA CORNACCHIA ET AL. Ecology,Vol. 99,No. 4 The average rainfall during the experiment was 1.12 mm/d, and there was no rainfall in 36 out of 49 d of the experiment (see Appendix S5). Growth rates were calculated in terms of shoot height as GR H = (H 2 ÀH 1 )/H 1 , with H 1 and H 2 being the shoot height (cm) on day 1 and day 49 of the experiment. Plant height and biomass are highly correlated for B. erecta and G. densa (e.g., Puijalon and Bornette (2004), Puijalon et al. (2005), Puijalon and Bornette (2006), and based on our previous sampling measurements in Appendix S6). Growth rate in height can be used as a nondestructive alternative to relative growth rate of biomass (Perez-Harguindeguy et al. 2013). Thus, we chose to assess plant size using plant height to minimize plant manipulation at the transplantation date. Moreover, this approach allowed us to keep transplantation time as brief as possible, which is important to avoid plant deterioration. In our case, the plants were harvested from within the same channel and immediately transplanted at the selected study locations without bringing them back to the laboratory for biomass measurements. Here, the initial transplanted individuals were referred to as "mother ramets". New ramets produced by mother ramets through vegetative reproduction were referred to as "daughter ramets", and stolons and daughter ramets together were defined as "juveniles". Shoot height, number of stolons, total stolon length, spacer length, and number of daughter ramets were measured on the transplants. Afterwards, biomass was separated into mother ramet and juveniles, dried in the oven at 60°for 48 h and weighed to obtain the dry mass of the transplants and the biomass investment in vegetative reproduction.
To characterize the flow velocity encountered by transplants for each treatment, both in the surroundings of C. platycarpa patches and on bare sediment, we measured flow velocities in the proximity of each transplant. Flow was measured for 100 s at 1 Hz using an Acoustic Doppler Velocimeter (ADV; FlowTracker, SonTek) at a water depth of 60% from the water surface, to obtain an estimate of average flow velocity over the water column. The study river was selected for its uniform channel structure (cross-section, water depth) and because it is artificially managed by the Compagnie Nationale du Rhône (CNR), maintaining stable conditions in terms of discharge and water levels all year round. Previous measurements at the study site showed that summer flow velocities were stable over time, and this trend was confirmed in the following summer (see Appendix S4: Fig. S1). Thus, flow velocity measurements were taken once during the experimental period to characterize the typical flow conditions in different locations around Callitriche platycarpa patches. The relative differences in velocity among treatments were assumed to be reasonably constant over time, despite some fluctuation in discharge. The flow velocities encountered by each transplant were subsequently correlated to their growth rates, survival and traits of vegetative reproduction at the end of the experiment.
One-way ANOVA was applied to test for significant differences in dry biomass of transplants between positions around existing patches. Post-hoc comparisons were performed using a Tukey HSD test. Survival of transplants between treatments was analyzed using Kaplan-Meier survival analyses and Mantel-Cox log rank tests with Bonferroni correction. The relationships between flow velocity and height increase, spacer length, daughter ramet dry mass, and between mother and daughter ramet height, were tested with a linear regression model. All statistical analyses were performed in R 3.1.2.

Model simulation of species coexistence patterns
Model simulations showing self-organized pattern formation demonstrated that scale-dependent feedbacks between macrophytes, sedimentation, and hydrodynamics could generate the patchy vegetation distribution observed in the field ( Fig. 2A). Regular patterns of vegetation, consisting of welldefined high biomass patches alternating with bare sediment with little vegetation, develop at intermediate flow velocities. The patches are streamlined and oriented in the main direction of the flow. Due to a scale-dependent interaction of vegetation with water flow, increased flow resistance locally reduces flow velocities within the vegetation, while water flow is diverted and accelerated between the vegetation patches (arrows in Fig. 2A). Sedimentation is promoted within the patches, up to a point where high sediment accumulation on the downstream side of the patches limits their further length growth in the streamwise direction. Our model highlights that self-organization processes between vegetation growth and hydrodynamics are a potential explanation for the patchy characteristics of many streams, especially at intermediate flow velocities.
When the pattern forming facilitator species P f is allowed to interact with the non-pattern forming beneficiary species P b , coexistence is promoted. A beneficiary species P b1 with low resistance to hydrodynamic stress is able to colonize the sheltered, low-flow areas in the wake region downstream of the P f patches, but is outcompeted within the patches themselves (Fig. 2B). A beneficiary species P b2 with lower growth rate r and higher resistance to hydrodynamic stress can coexist inside and locally around the margins of P f patches, near the high-flow areas created on the sides (Fig. 2C). Hence, our model shows that, in hydrodynamically stressful habitats, species with different resistance to flow stress can coexist through different spatial patterns, either in the wake of the patterned facilitator species P f , or locally inside and along the margins of the dominant patterns. These new niches are created by the hydrodynamic heterogeneity resulting from the self-organization process.
Our model analyses also highlight that the presence and strength of the interactions between facilitator and beneficiary species depend strongly on hydrodynamic conditions. The realized biomass of each species under homogeneous conditions (Eq. 5) shows that changes in incoming flow velocity determine the shift from dominance of one species, to stable coexistence, to dominance of another species (realized biomass distributions in light green, dark green and orange; Fig. 2D). At low incoming flow velocity (U in ), P b1 is the most successful competitor (Fig. 2D); as flow velocity increases, P b1 and P f can coexist within the range 0.07 ≤ U in ≤ 0.09. As incoming flow increases further, P f becomes the dominant species, until a range where it coexists with P b2 . At the highest flow velocities, P b2 is the most successful competitor due to its higher resistance to flow stress. Based   Vol. 99,No. 4 on the species realized niches along the flow velocity gradient, our model analysis also shows that in a spatial model for a given U in , a uniformly distributed P f would attenuate incoming flow velocity U in to a single realized velocity U e that would be more favourable for its growth. This flow velocity falls in the range where P f is predicted to be the only dominant species (Fig. 2D). Instead, in a spatial model for the same flow velocity, a self-organizing P f would separate the incoming flow into areas with low velocity (U e , inside the patches and in their wake) and areas with high velocity (U e , next to the patches), thus promoting coexistence and diversity by creating a much wider range of hydrodynamic conditions that provide the niches where each species can be dominant (Fig. 2D, E). Testing for hydrodynamic heterogeneity under self-organization highlights the very wide range of hydrological niches created by this process in the spatial model (Fig. 2E). The frequency distribution of flow velocities over the simulated domain shows that self-organization creates a much wider range of hydrodynamic conditions, compared to homogeneous environments. Self-organized patterning leads to a bimodal distribution of flow velocities, with a low-flow peak in vegetated areas, and a high-flow peak in unvegetated areas between plant patches (frequency distributions in dark green and blue; Fig. 2E). The self-organizing species therefore provides a spatial flow velocity gradient: low stress areas where less resistant species are more successful, and higher stress areas where more resistant species are dominant. Such hydrodynamic heterogeneity promotes coexistence by allowing all outcomes of species interactions to occur in space. Depending on the incoming flow velocity U in set at the beginning of the simulation, and on the species included in the model, the extent of the flow attenuation within the patches and acceleration around them (i.e., the ranges of realized velocity U e ) might be different (Fig. 2E). Our model highlights that, under self-organization, beneficiary species can persist in environments they would not normally inhabit based on average flow conditions. Therefore, facilitation expands the niches of the beneficiary species and allows them to withstand stronger hydrodynamic stress levels.

Comparison between simulated and observed species coexistence patterns
Spatial autocorrelation analysis to test for self-organization in the field shows that the spatial patterns of P f predicted by our numerical model display significant positive autocorrelation up to 1.5 -2 m distance, followed by significant negative autocorrelation at a distance up to 3 - There is a clearly observable agreement between the spatial correlation function from the field patterns of C. platycarpa and the results of the autocorrelation analysis on the predicted patterns. Obviously, differences in patch geometry between the model and the real-world patches appear upon visual inspection (Fig. 4A, B), as the model only captures a subset of the relevant processes. Yet, the spatial analysis reveals the regularity of the spatial pattern, with plant aggregation on short scales (positive autocorrelation) and overdispersion (negative autocorrelation) at larger scales. The mean wavelength of the spatial patterns is, however, different: C. platycarpa patches are located every 5 m in the model and 8 m in the field. Autocorrelation analysis of C. platycarpa patches from our aerial pictures either showed significant positive autocorrelation up to 2 m distance, followed by significant negative autocorrelation from 3 to 5 m ( Fig. 4B; black line in Fig. 4D), or it showed a directional effect of significant positive autocorrelation up to 6 m distance, but without negative correlation at any distance due to merging of neighbouring patches ( Fig. 3B; black line in Fig. 3D). Hence, in the first case (Fig. 3D) we found streamlined bands of vegetation distributed in the direction parallel to the main flow direction, with no clear gap between the patches due to their merging. In the second case (Fig. 4D), we found regular vegetation patches oriented parallel to the main flow direction, at a distance of roughly 8 m from each other.
When a second species P b is included in our model, the predicted outcome of species interaction is that P b can coexist in the low-flow areas created in the wake of the patches of the pattern-forming species P f (Fig. 3A). Spatial cross-correlation analysis of P f with P b indeed shows a significant positive association of the beneficiary species in the wake of existing patches of the facilitator, as shown by the positive peak in the cross-correlation coefficient at around 1.0 m distance from and accelerated outside (indicated by arrow size and color, from yellow to red). (B) Beneficiary species characterized by low resistance to hydrodynamic stress (light green) colonize the sheltered, low-flow areas in the wake of the P f patches (dark green), while being outcompeted within the patches themselves. (C) Beneficiary species with lower growth rate and higher resistance to hydrodynamic stress (orange) can coexist inside and locally around the P f patches (dark green), near the high-flow channels created next to them. (D) Realized niches of P f , P b1 and P b2 along the hydrodynamic stress gradient in the homogeneous model. Dashed lines indicate the limits between the flow velocity ranges where either one species is dominant, or two species coexist. In a spatial model for a given U in , a uniformly distributed P f would attenuate incoming flow velocity U in to a single realized velocity U e . This flow velocity falls in the range where P f is predicted to be the only dominant species (based on the species realized niches along the flow velocity gradient, in D). Instead, for the same flow velocity, a self-organizing P f would separate the incoming flow into areas with low velocity (U e , inside and downstream of the patches) and areas with high velocity (U e , next to the patches), thus creating a wider range of hydrodynamic conditions that provide the niches where each species can be dominant (in D, E). Parameters used are r f = 1.19, a fb1 = 0.6, a b1f = 1.42, k b1 = 390, r b1 = 0.94, a b2f = 0.83, k b2 = 100. Other parameters as in Table 1. (E) Hydrodynamic heterogeneity generated by self-organization in the spatial model: frequency distribution of depth-averaged flow velocities within vegetated (dark green) and unvegetated cells (blue) of the simulated domain. The two subfigures refer for the two beneficiary species: P b1 (top figure) and P b2 (bottom figure).
( Fig. 2. Continued) April 2018 SELF-ORGANIZED LANDSCAPE OF FACILITATION 839 them (blue line in Fig. 3C). Parallel to the main flow direction, this spatial cross-correlation function shows correspondence to the species coexistence patterns found in the field. Berula erecta showed a significant positive association in the wake of C. platycarpa patches (Fig. 3D). Our analysis shows the peak of the beneficiary species around the same downstream location in the field (1.5 m; Fig. 3D) and simulations (1 m; Fig. 3C). In contrast, the cross-correlation analysis in the direction perpendicular to the main flow reveals a difference in behaviour between the simulated and observed Individual patches can be obscured because they can merge and grow above another, but 7 patches of Berula and more than 20 of Callitriche were present in the reach. Please note the scale difference compared to the model in (A). Auto-and cross-correlation functions of species distribution patterns from model simulations (C) and field observations (D) in the direction parallel to the main water flow. Auto-and cross-correlation functions of species distribution patterns from model simulations (E) and field observations (F) in the direction perpendicular to the main water flow. In C and E, black lines are the autocorrelation functions for the simulated spatial patterns of P f ; blue lines are the cross-correlation functions between P f and P b . In D and F, black lines are the autocorrelation functions for Callitriche platycarpa; blue lines are the cross-correlation functions between Callitriche and Berula. Closed dots represent significant values.

840
LORETA CORNACCHIA ET AL. Ecology,Vol. 99,No. 4 patterns. Field observations show that Berula is located along the outer edges of the Callitriche patches (blue line in Fig. 3F). This pattern differs from the model results, where the beneficiary species occupies the region immediately downstream of the facilitator patches (Fig. 3E).
When P b is used to model a species with higher resistance to flow stress, a different pattern of coexistence is observed: the beneficiary species grows both within the patches and in the open interspaces around the pattern-forming species ( Fig. 4A; blue line in Fig. 4C). This predicted pattern of Individual patches can be obscured because they can merge and grow above another, but 9 patches of Groenlandia and more than 30 of Callitriche were present in the reach. Please note the scale difference compared to the model in (A). Auto-and cross-correlation functions of species distribution patterns from model simulations (C) and field observations (D) in the direction parallel to the main water flow. Auto-and cross-correlation functions of species distribution patterns from model simulations (E) and field observations (F) in the direction perpendicular to the main water flow. In C and E, black lines are the autocorrelation functions for the simulated spatial patterns of P f ; blue lines are the cross-correlation functions between P f and P b . In D and F, black lines are the autocorrelation functions for Callitriche platycarpa; blue lines are the cross-correlation functions between Callitriche and Groenlandia. Closed dots represent significant values.

April 2018
SELF-ORGANIZED LANDSCAPE OF FACILITATION 841 coexistence is in strong agreement with field observations on coexistence patterns of Groenlandia densa and Callitriche platycarpa, where Groenlandia tended to coexist within and along the margins of Callitriche patches ( Fig. 4B; blue line in Fig. 4D). In both cases, the two species are positively associated up to 2 m distance (i.e., where the patches of the patterned species are located), but negatively or non-significantly correlated from 2 to 5 m distance (i.e., where the patterned species is absent due to the negative feedback on its growth). The relationship between Callitriche and Groenlandia in the direction perpendicular to the flow in the field still shows a pattern of coexistence (Fig. 4F), as confirmed by the analysis of the model predictions (Fig. 4E), while also highlighting a shift in the lateral distribution of the two species as Groenlandia tends to grow along the margins of Callitriche patches. showed a significantly higher increase in shoot height compared with transplants on the "Bare sediment" treatment (t-test, t = 4.3, df = 4.387, P = 0.02 for Berula; t = 5.5, df = 1.839, P = 0.04 for Groenlandia). The intensity of this effect was correlated with the reduction in flow velocity created by the facilitator species (r 2 = 0.96, P = 0.0004 for Berula, r 2 = 0.82, P = 0.03 for Groenlandia; Fig. 5).
Vegetative reproduction.-No difference in dry mass invested in vegetative reproduction was found for either species between transplant positions (Fig. 6C, F). Dry mass investment was not correlated with incoming flow velocity for B. erecta (r 2 = 0.0168, P > 0.05) or for G. densa (r 2 = 0.48, P = 0.19). A significant negative correlation was found between the average spacer length in the transplants and incoming flow velocity for B. erecta (r 2 = 0.84, P = 0.01; Fig. 6A). The correlation was not significant for G. densa (r 2 = 0.64, P = 0.19; Fig. 6D). A significant positive correlation was found between the height of the mother ramet transplant and their average daughter ramet height for both B. erecta (r 2 = 0.65, P = 0.05) and G. densa (r 2 = 0.85, P = 0.02) (Fig. 6B, E).

DISCUSSION
In a combined mathematical and empirical study, we reveal that bio-physical feedbacks between in-stream submerged plants and streamflow can generate spatial heterogeneity in hydrodynamic conditions that create new niches, promoting species coexistence in streams. Central to this landscape of facilitation is spatial self-organization of submerged aquatic vegetation by means of deflection of water flow by the facilitator species, Callitriche platycarpa, which generates a patterned landscape of Callitriche patches. Our mathematical model shows that (1) the hydrodynamic heterogeneity results from the self-organization process and (2) it promotes coexistence by creating new niches for species that are adapted to a wider variety of environmental conditions. Species distribution patterns from our numerical model showed similarities with the spatial aggregation of macrophyte species around Callitriche platycarpa patches observed in the field at the reach scale. A field transplantation experiment revealed that species coexistence results from a positive interaction due to stress amelioration, as the growth of these beneficiary species was facilitated by the FIG. 6. Relationships between flow velocity within and around Callitriche platycarpa patches, and traits of vegetative reproduction for Berula erecta (A, C) and Groenlandia densa (D, F) at the end of the experiment (t = 49 d). Relationship between mother and daughter ramet height for Berula erecta (B) and Groenlandia densa (E). hydrodynamic stress reduction mediated by Callitriche patches. Moreover, the effects of self-organized pattern formation on species interactions go beyond the spatial structure of the vegetation community. By affecting clonal growth traits, Callitriche patches also affect the density of the patches of other species, and therefore the spatial organization and appearance of vegetation patterns for the beneficiary species. Our study highlights that species coexistence in streams is, in part, explained by a biophysical feedback process that creates a heterogeneous landscape offering facilitative effects.

Landscapes of facilitation through self-organized patchiness
Current theory largely ignores the spatial dimension when considering facilitative effects between species (Callaway (2007), Smit et al. (2007), Cavieres et al. (2014); but see van de Koppel et al. (2006van de Koppel et al. ( , 2015 for a review). Facilitative interactions are for the most part considered within the tussocks or patches of the facilitator species, and to date experiments have focused on this local scale, as beneficiary species are mainly considered to be living inside the facilitator patches (e.g., nurse plants in drylands; Callaway and Walker (1997), Badano and Cavieres (2006); but see Pescador et al. (2014)). Through this approach, many studies have shown the importance of facilitation but few have looked at its spatial variability. Here, we reveal that in self-organized ecosystems, facilitative interactions are far from being homogeneous in space, and display strong spatial heterogeneity due to the balance between positive and negative feedbacks. The selforganizing process leads to spatial separation of competition and facilitation, with opposite effects balancing throughout the landscape. Similar long-distance effects through modification of physical forcing by ecosystem engineers have also been observed in other systems, such as mussel beds on tidal flats (Donadi et al. 2013) or between adjacent tropical ecosystems at the landscape scale (Gillis et al. 2014). The heterogeneity of facilitation and its spatial effects are important processes that have been identified in previous studies (Bruno 2000, Bruno and Kennedy 2000, van de Koppel et al. 2006, although not in the context of self-organized ecosystems. Hence, we show that self-organization acts as a strong structuring force of community composition and distribution by creating spatial variability in environmental conditions, leading to facilitative interactions at different spatial scales.
Our results emphasize that by triggering a self-organized pattern, a single engineering species may create a "landscape of facilitation", where multiple mechanisms of coexistence co-occur due to the conditions created by the self-organized process. The conditions include: low stresshigh competition inside the patch; low stresslow competition downstream of the patch; and high stresslow competition next to the patch. As the facilitative effects described here extend over longer distances, species with higher resistance to stress can locally colonize the open interspaces around the patches, exploiting the new niches created by the negative feedback without being exposed to high competition; less tolerant species can grow at a certain distance from the patch, where the positive feedback of stress reduction is still present, but there is no negative effect of competition.
The comparisons between field vegetation patterns and model outputs highlight areas for further detailed experiments and model improvement. Our model is minimalistic and does not capture all of the relevant processes that occur in real streams. For instance, as the vegetation density increases, canopy-scale turbulence can lead to higher sediment resuspension within the vegetation (Yang et al. 2016), creating patterns of enhanced or diminished turbulence and sediment deposition in different locations. Moreover, the scaling of stem-scale and patch-scale turbulent wakes can limit the deposition of fine material downstream of a patch (Chen et al. 2012, Liu and. These findings suggest that Berula erecta might occupy optimal zones where sediment can be deposited, downstream of the termination of the turbulent wake structure and along the outer edges of Callitriche patches. Consistent with this, earlier studies by Sand-Jensen (1998) observed turbulent eddies and sediment erosion at the rear end of macrophyte patches of species with an overhanging canopy. These complex patterns in turbulence and sediment deposition are interesting possible extensions of the model that will provide an even more elaborate mechanistic basis for habitat and species diversity in streams.
Although our model depicts a simplification of the complex hydrodynamic-vegetation interactions, the comparison between the predicted and observed spatial patterns suggests that the spatial distribution of Berula erecta is similar to that of a beneficiary species with lower resistance to hydrodynamic stress, while Groenlandia densa exhibits greater behavioural similarity to species with higher resistance to stress. The differences in stress resistance between the two species are also supported by our transplantation experiments. For Groenlandia densa, we found a steeper slope and larger yintercept of the negative relationship between flow velocity and growth rate, compared to Berula erecta (Fig. 5). As the regression line for Groenlandia is located above the line for Berula across the whole range of flow velocities in our experiment, the former appears to perform consistently better in response to flow stress. Survival results for Berula erecta showed significantly higher mortality within the patch than in the other treatments, suggesting that short-range competition for light prevails in that location. However, while we found a facilitative effect in terms of growth rates of the initial transplanted individuals, we found no effect on the biomass they invested in vegetative reproduction (through clonal growth). This observation is consistent with the ability of B. erecta to maintain its investment in vegetative growth and produce a more compact clonal growth form, despite the increased flow stress (Puijalon et al. 2005, Puijalon andBornette 2006). Therefore, self-organization processes allow the coexistence of species with a wide range of growth strategies and sensitivity to stress.

Effects of self-organization on species coexistence
The process of pattern formation allows species to coexist, even if the number of resources on which they grow would predict competitive exclusion (Gause 1932). The results from our study on submerged macrophytes in streams are in accordance with the only known previous theoretical studies of pattern formation and species coexistence albeit on arid savannas (Gilad et al. 2004, Baudena andRietkerk 2013, 844 LORETA CORNACCHIA ET AL. Ecology, Vol. 99, No. 4 Nathan et al. 2013). However, while these studies found coexistence of two species within the same spatial pattern (i.e., overlapping patches), we found that self-organization effects act both locally and at distance beyond the limits of the facilitator canopy (in the order of a few meters of the river reach in our study). Hence, self-organization can provide a potential explanation for the high biodiversity observed in many natural communities, despite theoretical predictions of low species coexistence. Self-organization differs from ecosystem engineering and local facilitation between species in important ways. While ecosystem engineering creates a local positive feedback, selforganized patchiness also results from both a positive and a strong negative feedback. This negative feedback has a twofold role. First, it prevents the facilitating species from dominating the entire habitat. Second, it changes environmental conditions within the inter-patch spaces, allowing for the coexistence of a wide range of species as compared to the original, more homogeneous habitat. Therefore, the emergence of self-organized patterns produces distinct spatial signatures in plant community structure that might be discerned from local facilitation effects.
The creation of new niches and the effects on biodiversity arising from facilitation can benefit both plant and animal species. For instance, fish can use both the shelter provided by plants as protection from predation, and the high-flow areas around patches as spawning and feeding grounds (Kozarek et al. 2010, Marjoribanks et al. 2016; and suspension-feeding invertebrates (e.g., blackfly larvae) can grow on the edge of submerged macrophyte patches, such as Ranunculus sp. where higher current velocities increase the flux of resources (Wharton et al. 2006). Thus, spatial self-organization has the ability to affect many species within stream communities at different trophic levels.

Relevance beyond stream ecosystems
The importance of pattern formation in promoting species coexistence is likely to be relevant for a wide range of selforganized ecosystems. In many of these systems, at least one habitat-forming species provides structure for an entire community. For example, periodic vegetation patterns in arid or semi-arid systems create different levels of edaphic and climatic stress for other species (Couteron 2001, Rietkerk et al. 2002. In coastal environments, mussel beds on relatively homogeneous intertidal flats reduce wave stress and increase habitat structural complexity and species richness (Guti errez et al. 2003, van de Koppel et al. 2005, Donadi et al. 2013, Christianen et al. 2016) and salt marsh plants create different spatial patterns of sediment deposition, salinity and redox conditions (Howes et al. 1980, Callaway 1994, Hacker and Bertness 1999. Thus, as self-organized patterns emerge as a widespread phenomenon, landscapes of facilitation which enhance species coexistence and biodiversity are likely to be of similar ecological importance.
In ecosystems with limited underlying heterogeneity in abiotic conditions, self-organization acts as a powerful structuring force of community composition and distribution. These findings can be used to inform ecological restoration projects, which aim to maximize biodiversity through the preservation or re-introduction of self-organized species. Exploring the implications of species coexistence promoted by self-organization on food web structure is also an interesting topic for future studies. Understanding of the intricate way in which competition and facilitation interact in many ecosystems is key to successful management of their biodiversity.