Similar fish species composition despite larger environmental heterogeneity during severe hypoxia in a coastal ecosystem

Abstract Environmental heterogeneity is one of the most influential factors that create compositional variation among local communities. Greater compositional variation is expected when an environmental gradient encompasses the most severe conditions where species sorting is more likely to operate. However, evidence for stronger species sorting at severer environment has typically been obtained for less mobile organisms and tests are scarce for those with higher dispersal ability that allows individuals to sensitively respond to environmental stress. Here, with the dynamics of fish communities in a Japanese bay revealed by environmental DNA metabarcoding analyses as a model case, we tested the hypothesis that larger environmental heterogeneity caused by severe seasonal hypoxia (lower concentration of oxygen in bottom waters in summer) leads to larger variation of species composition among communities. During summer, fish species richness was lower in the bottom layer, suggesting the severity of the hypoxic bottom water. In contrast to the prediction, we found that although the environmental parameters of bottom and surface water was clearly distinct in summer, fish species composition was more similar between the two layers. Our null model analysis suggested that the higher compositional similarity during hypoxia season was not a result of the sampling effect reflecting differences in the alpha or gamma diversity. Furthermore, a shift in the species occurrence from bottom to surface layers was observed during hypoxia season, which was consistent across species, suggesting that the severe condition in the bottom adversely affected fish species irrespective of their identity. These results suggest that larger environmental heterogeneity does not necessarily lead to higher compositional variation once the environmental gradient encompasses extremely severe conditions. This is most likely because individual organisms actively avoided the severity quasi‐neutrally, which induced mass effect‐like dispersal and lead to the mixing of species composition across habitats. By showing counter evidence against the prevailing view, we provide novel insights into how species sorting by environment acts in heterogeneous and severe conditions.

geneity caused by severe seasonal hypoxia (lower concentration of oxygen in bottom waters in summer) leads to larger variation of species composition among communities. During summer, fish species richness was lower in the bottom layer, suggesting the severity of the hypoxic bottom water. In contrast to the prediction, we found that although the environmental parameters of bottom and surface water was clearly distinct in summer, fish species composition was more similar between the two layers. Our null model analysis suggested that the higher compositional similarity during hypoxia season was not a result of the sampling effect reflecting differences in the alpha or gamma diversity. Furthermore, a shift in the species occurrence from bottom to surface layers was observed during hypoxia season, which was consistent across species, suggesting that the severe condition in the bottom adversely affected fish species irrespective of their identity. These results suggest that larger environmental heterogeneity does not necessarily lead to higher compositional variation once the environmental gradient encompasses extremely severe conditions. This is most likely because individual organisms actively avoided the severity quasi-neutrally, which induced mass effect-like dispersal and lead to the mixing of species composition across habitats. By showing counter evidence against the prevailing view, we provide novel insights into how species sorting by environment acts in heterogeneous and severe conditions.

| INTRODUC TI ON
Understanding the processes controlling the variation in species composition across multiple habitats has been a central challenge in ecology. The metacommunity framework (Leibold et al., 2004;Thompson et al., 2020) provides a way to integrate multiple community assembly processes and roles played by environmental heterogeneity and dispersal. The species sorting paradigm, one of the paradigms in the metacommunity framework, predicts that species composition deterministically aligns with environmental gradients as determined by each species' requirements for some abiotic conditions (Grinnell, 1917;Tilman, 1982). In this paradigm, species composition is expected to be more dissimilar if environmental heterogeneity among habitats is larger (MacArthur et al., 1966;Ricklefs, 1977). Meanwhile, the mass effect paradigm considers the higher rate of dispersal of individuals, which promotes the mixing of species composition among localities (Mouquet & Loreau, 2002). In this case, similar species composition may be observed even among environmentally heterogeneous habitats. The neutral dynamics paradigm, where species are modeled as equivalent in their fitness and competitive ability, poses a prediction similar to the mass effect that the environment-species composition relationship is no longer detectable. These paradigms are not discretely separated; rather, they are considered as the continuum determined by the dispersal rate and the degree of interspecific differences (Thompson et al., 2020).
In this context, much effort has been devoted to understanding the condition in which the relative importance of species sorting by environment over the mass effect or neutrality varies Gravel et al., 2006).
Growing body of literature suggests that species sorting plays a more important role when the abiotic environment is more severe (Chase, 2007;Daniel et al., 2019;Guo et al., 2014;Lepori & Malmqvist, 2009). This is because at the severe environment, where the growth rate of species is lowered or the mortality rate is high, only a limited number of species that possess traits adapted to the severe environment can survive. This leads to a prediction that compositional variation becomes great when environmental heterogeneity is large enough to encompass the severe environment where species sorting is more likely to act ( Figure 1a). However, there are two major knowledge gaps in the previous literature, which hinders our ability to fully understand the relationship between environmental severity and community assembly.
First, the current understanding assumes that responses to environmental severity differ among species. However, a limited number of studies indicated exceptional cases in which species difference becomes small and species sorting is weak at the severest conditions (Kim, 2019;Lepori & Malmqvist, 2009). For example, Lepori and Malmqvist (2009) found that although moderate disturbance in streams (frequent flooding) filtered out some species that are intolerant of the imposed stress, too severe disturbance did not work as species sorting because it induced local extinction of individuals irrespective of species traits. Therefore, K E Y W O R D S community assembly, eDNA, environmental filtering, environmental heterogeneity, fish, hypoxia

T A X O N O M Y C L A S S I F I C A T I O N
Community ecology F I G U R E 1 Two alternative hypotheses regarding the effect of environmental heterogeneity on variation in species composition. (a) Habitats with distinct environments (e.g., different water layers) have more dissimilar species composition when the environmental heterogeneity is large because of species sorting. (b) Even habitats with distinct environment share similar species composition because the severest conditions (e.g., extreme hypoxia at the bottom) induce higher mortality rates or movement toward preferable habitats (e.g., the surface layer) irrespective of species identity the effect of environmental severity on the strength of species sorting depends on whether the species difference is maintained at the most stressful conditions, but little attention has been paid for this possibility in the current framework (Lepori & Malmqvist, 2009). The second knowledge gap is that the evidence that species sorting is more likely at the severe environment typically comes from systems where dispersal of individuals is limited or passive, such as plants or aquatic organisms in ponds (e.g., Chase, 2007;Guo et al., 2014), and little is known for organisms with higher dispersal ability. Unlike less mobile organisms, highly mobile organisms can actively escape from unpreferable severe environments, which can lead to the concentration of many species in preferable habitats, regardless of their ability to tolerate the stress. In that case, while the environmental severity indeed affects organisms' distribution, it induces the mass effect-like dispersal and no longer operates as 'species' sorting. Therefore, when the environmental heterogeneity becomes large enough to encompass the severest conditions, it is probable that the compositional variation would not be larger or even be smaller despite the environmental heterogeneity ( Figure 1b), due to smaller species differences (i.e., more neutrality) or nonspecies-specific movement toward preferable habitats (i.e., mass effect-like dispersal).
In marine ecosystems, one of the environmental conditions that most drastically influence community structures is hypoxia or the decline in dissolved oxygen (DO; Breitburg et al., 2018;Diaz & Rosenberg, 1995, 2008. One of the typical forms of hypoxia is seasonal: in summer, water warming forms the stratification of water columns, which leads to the environmental heterogeneity along water depth. Specifically, the reduced concentration of DO in the bottom water to a critical level imposes extreme stress for demersal organisms, such as increased mortality rates (Hrycik et al., 2017;Levin, 2003). Therefore, the species-sorting paradigm ( Figure 1a) expects that species that do not possess morphological or physiological traits adapted to the hypoxic conditions are likely to be filtered out from the bottom water (Bickler & Buck, 2007;Kodama & Horiguchi, 2011;Levin, 2003), and different species composition between bottom and surface water is expected during that season. In contrast, it is also probable that the bottom hypoxia was so intense that many species suffer from it or actively avoid it to a similar extent (Eby & Crowder, 2002;Pihl et al., 1991), and species can be apparently neutral with respect to the responses to the severe hypoxia, weakening the effect of species sorting ( Figure 1b).
The objective of this study is to explore whether among- (1) species composition is more dissimilar between bottom and surface layers during the hypoxia season (summer), because the environmental heterogeneity is larger, and the bottom environment is extremely stressful which strongly sorts species that are adapted to the severity ( Figure 1a). Alternatively, (2) compositional difference between bottom and surface waters is not evident or even smaller during hypoxia season despite larger environmental heterogeneity. This is because the extremely severe environment in the bottom negatively affects fish occurrences irrespective of the species identity and thus the severity does not operate as 'species' sorting ( Figure 1b), which leads to the mixing of species composition. To this end, we took advantage of the availability of data of seasonal dynamics of fish communities in Tokyo Bay revealed by the environmental DNA (eDNA) metabarcoding analysis (Hongo et al., 2021) which has emerged as a powerful and efficient tool for monitoring biodiversity especially in aquatic ecosystems (Bohmann et al., 2014;Deiner et al., 2017). By extracting genetic materials shed from organisms from environmental samples (e.g., water and soils), this technique allows for the noninvasive assessment of biodiversity, with a detection ability comparable or superior to that of traditional sampling methods (Sigsgaard et al., 2017;Thomsen et al., 2012;Yamamoto et al., 2017). Importantly, owing to the short persistence (from a few hours  to several days [Thomsen et al., 2012]) and low diffusion rate (less than 100 m;

| Environmental variables
Three environmental variables, namely concentration of DO, water temperature, and salinity, were measured at the same time as the water sampling. During summer, the water temperature and DO concentration were lower in the bottom than in the surface (Figures 2a and S1b). Specifically, from June to September, the minimum value of DO concentration was lower than 2 mg/L ( Figure   S1b), a common threshold of hypoxia in previous studies (Breitburg et al., 2018;Hrycik et al., 2017;Pihl et al., 1991). Following the established threshold in the literature, we determined these months as "hypoxia" season and other months as "normoxia" (normal levels of oxygen). Note that even in the hypoxia season, the DO values considerably varied among samples, with many being higher than the threshold ( Figure S1b). Nevertheless, we classified all samples of the months when the minimum DO value falls below the threshold as the hypoxia season. This is because analyzing the samples obtained in the same period as one unit allows us to evaluate the effects of environmental heterogeneity on the compositional variation among samples.
In addition, while significant decreases in DO concentration in summer were observed at most of the sites, this was not the case at sites near the mouth of the bay ( Figure S1) because of deep water depth and water exchange with the ocean (Fujiwara & Yamada, 2002), which may unexpectedly blur the results. Nevertheless, our main finding was unchanged when we excluded the sampling sites near the bay mouth from the analyses ( Figure S2).

| Data analysis
To minimize unexpected biases, we filtered out the species that did not meet the following two criteria from further analyses: (1) those whose names were not matched during the search in Fishbase, a global database of fish (Froese & Pauly, 2019; https://www.fishb ase. se/), and (2) those whose habitats were known to be out of Tokyo Bay based on information from a range of studies (Kouno et al., 2011;Nakabo, 2002). The first criterion mainly filtered out hybrid species.
The second one, which was based on each fish species' distribution in Japan and observation records in Tokyo Bay, removed fish species known to be freshwater or endemic in other regions which were detected likely due to sequencing errors or contamination from rivers or nearby fish markets.
All statistical analyses were performed using R version 4.0.2 (R Core Team, 2020) and the figures were illustrated using the ggplot2 package (Wickham, 2016). and 41.1% of the total variation respectively, and the first two MDS axes of the species composition. Differences in the environment and composition among seasons and layers were tested using permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001), based on Euclidean and Jaccard distance metrices for the environment and composition data, respectively. We included the interaction term between season and layer, as well as their independent effects, in the PERMANOVAs to determine whether there exists a seasonal difference in the extent of surface-bottom divergence. Since the number of sampling sites was small (n = 2) due to logistical problems, we did not conduct the above analyses for the October data. The PERMANOVAs were implemented by the adonis function in the vegan package (Oksanen et al., 2019).
Results of the PERMANOVA that tested the interactive effect of season and layer on species composition may be biased by seasonal differences in the size of regional species pool (gamma-diversity) and local species richness (alpha-diversity), because higher gammadiversity and lower alpha-diversity generally lead to larger compositional difference simply due to the sampling effect (Kraft et al., 2011;Mori et al., 2015). To examine this undesired possibility, we relied on a null model analysis. Our null model randomly shuffled the occurrence of species in each sample with total number of occurrences of each species and the samples' species richness fixed.
Note that the randomization was conducted separately for each month as a unit sharing the same regional species pool. For each month, 999 randomized community matrices (sample × species) were generated based on the matrix-swap algorithm (Gotelli, 2000).
The randomization was implemented by the RandomizeMatrix in the package picante (Kembel et al., 2010). For each pair of samples, To examine the species-specific responses to bottom hypoxia, we compared the frequency of occurrence of fish species between seasons (hypoxia or normoxia) and water layers (surface or bottom), using a GLMM with species-specific random intercept and coefficients. The response variable was the presence probability of species i in a sample s (p i,s ), and the explanatory variables were season, layer and their interaction: The season and layer were transformed into dummy variables (hypoxia: 0, normoxia: 1, and surface: 0, bottom: 1). To assure sufficient statistical power, we restricted this analysis to the frequently observed species, that is, top 10% most frequently observed species

| RE SULTS
We obtained 308 water samples for the eDNA extraction, and gene amplification were successful for 281 samples, from which about 28 M reads were obtained and approximately 4 M reads were assigned as fish species (Hongo et al., 2021). The gene amplifications were not detected from negative control filters. Among these, corresponding environmental data were available for 226 samples (n = 23, 22, 22, 21, 20, 20, 18, 21, 20, 2, 15, and 22 from January to December, respectively). A total of 220 fish species were detected in the 226 samples. Of these, our first criterion filtered out 12 species, and the second criterion eliminated 40 species, resulting in 168 species retained in the subsequent analyses (Table S1).
During normoxia, species richness per sample was higher in the bottom layer (mean: 10.1, standard deviation [SD]: 4.78, Figure S3) than in the surface layer (mean: 6.24, SD: 3.05). However, such a difference was small during hypoxia: mean species richness was 9.15   (Figures 4b and S4). The GLMM, where the occurrence of each fish species was modeled as a function of water layer, season, and their interaction with species-specific random intercept and slopes, showed that layer and season did not influ- Indeed, the estimated species-specific coefficients of the interaction term were positive for all the 17 species ( Figure S5). Furthermore, the likelihood test that compared the above GLMM with another model without the species-specific random effect for the interaction term (Table 1b), showed no significant differences between the models (χ 2 = 2.56, df = 1, p = .109). These results collectively imply that the negative effect of hypoxia at bottom on the occurrence probability can be regarded as consistent across species.
We found a significant interactive effect of water layer and water depth of sampling sites on species composition (PERMANOVA,  Figure S6). This suggests that the compositional difference between bottom and surface layers was smaller at shallower sites.

| DISCUSS ION
Species sorting, one of the metacommunity paradigms, predicts that environmental heterogeneity is a major source of compositional variation. Based on pervasive understanding that environmental filtering works more in abiotically severer habitats (Chase, 2007;Daniel et al., 2019;Guo et al., 2014;Lepori & Malmqvist, 2009), among-habitats compositional variation is expected to be greater  Note: (a) Species-specific random effects were modeled for the intercept, and slopes of the three explanatory variables. (b) Species-specific random effects were modeled for the intercept and slopes of independent effect of layer and season only (not for the interaction). Estimated coefficients, standard errors (SE), p-values of the fixed effects, estimated variance of the random effect are shown. Significant (p < .05) estimates are shown with bold. The likelihood ratio test indicates that the fits of these models were not significantly different (χ 2 = 2.56, df = 1, p = .109).

TA B L E 1
Results of the two GLMMs that modeled the presence of species as a function of water layer (surface or bottom), season (hypoxia or normoxia), and their interaction with different modeling for the random effects when the environmental heterogeneity is so large that the gradient encompasses the most severe habitats (Figure 1a). In the present study system, extremely low DO concentration (<2 mg/L) in the bottom water during summer ( Figure S1) served as the severe condition, decreasing the species richness in the bottom layer ( Figure S3). As tolerating such extreme oxygen conditions requires certain morphological and physiological traits (Bickler & Buck, 2007;Mandic et al., 2009), different sets of species were likely at surface and bottom layers. However, against this expectation, we found that species composition was more similar during the hypoxic season ( Figure 2).
The null model analysis indicated that this pattern was robust even after accounting for the effect of sampling bias from the species pool ( Figure 3).

| Why less compositional variation during hypoxia?
The result that species composition was more similar during hypoxia  et al., 1991). Likewise, Eby and Crowder (2002) found that while the threshold for avoiding hypoxia varied among 10 fish species, it was larger than 2 mg/L for all species. Therefore, it is likely that, although the decrease in DO in the bottom does serve as a filter for species composition, once it exceeds a certain level of severity, fish behaviors (e.g., survival, foraging and growth) will be adversely affected regardless of species identity, at which point limited oxygen no longer plays a role in species sorting. Indeed, in our results, most of the common fish species (14 out of 17) showed a shift of their occurrence frequency from the bottom to surface layer ( Figure 4).
Moreover, our GLMM indicated that the random variation of the coefficient of the interaction term between water layer and season was smaller than that of their independent effects (Table 1a), and the species-specific coefficients were in the same direction for all the common species ( Figure S4). In addition, the two GLMMs, with and without the random variation in the coefficient of the interaction term, yielded similar fittings. These results collectively imply that even though species differ in their occurrence probability at different layers and seasons (as shown by the estimated large variance of random effects, Table 1), responses to the bottom hypoxia (i.e., less frequent occurrence in the bottom layer in hypoxia) do not vary much across species. Earlier experimental and observational studies have suggested that, for some of the common species in our dataset, behavioral, physiological, or distributional responses to hypoxia were seen at the DO value similar to our threshold (2 mg/L): Engraulis japonicus (1.12-2.36 mg/L, Oda et al., 2018), Mugil cephalus (2 mg/L, Wannamaker & Rice, 2000), Sardinops sagax (2.0-3.0 mg/L, Kreiner et al., 2009), Scomberomorus niphonius (4.0 mg/L, Shoji et al., 2005), Trachurus japonicus (2 mg/L, Yamamoto, 1991), while information was not available for the other species. Partially supporting our conclusion that the effect of bottom hypoxia is nonspecies specific, a meta-analysis revealed that effects of hypoxia on 29 fish species were not well explained by species ecological characteristics (Hrycik et al., 2017). Our results and other literature collectively support our second hypothesis that species sorting is weak during hypoxia because extremely low DO concentration harms fish irrespective of the species identity.
The observed shifts in the relative occurrence frequency from the bottom to surface waters during hypoxia can be a result of upward movements of individuals or higher mortality in the bottom (Diaz & Rosenberg, 2008;Levin, 2003). As fish are highly mobile, they can respond sharply to the environmental changes by avoiding unfavorable conditions. Therefore, when the bottom water was hypoxic in summer, they could move upward to escape from the stressful habitat (Eby & Crowder, 2002;Pihl et al., 1991). This would be particularly true for small pelagic fish such as those frequently observed in this study. Moreover, we observed increases in species richness per sample in the surface water during hypoxia, as well as the decreases in the bottom ( Figure S3), suggesting that species at the bottom layer moved up to the surface during hypoxia.
Meanwhile, evidence seems relatively weak for the case that higher mortality in the bottom is responsible for the observed shifts in the species occurrences. If so, species richness should have only decreased in the bottom and not increased in the surface, contrary to our result ( Figure S3). Nevertheless, some species, especially demersal fish, may have experienced higher mortality rate. For example, Branchiostegus japonicus whose habitat is specialized to the bottom (M. Ishii & H. Mita, personal communication) almost disappeared during the hypoxia period ( Figure S4) and its occurrence did not recover after the hypoxic condition improved, suggesting its higher mortality during hypoxia. Although being beyond the scope of this paper, comparing the responses to bottom hypoxia between demersal and small pelagic fish species would be a fruitful direction of future studies.

| Implications for community assembly studies
Our results have general implication in the context of community assembly. We found that fish compositional difference between surface and bottom layers in a bay was smaller during hypoxia season when the environmental heterogeneity was large. This was likely because the extremely hypoxic condition at the bottom induced masseffect like distributional shifts from the surface and surface layers irrespective of the species identity, or fish are similarly adversely affected by the severity and thus species differences are small (i.e., more neutral) at the bottom. Although the result appears as conflicts with the prevailing view that environmental heterogeneity leads to larger compositional variation and species sorting plays a more important role in stressful habitats (Chase, 2007;Daniel et al., 2019;Guo et al., 2014;Lepori & Malmqvist, 2009), a limited number of studies agrees with our finding. Lepori and Malmqvist (2009)

| Caveats and utility of eDNA metabarcoding analyses
Although the ability of eDNA metabarcoding analyses to detect species from environmental samples has been established (Sigsgaard et al., 2017;Thomsen et al., 2012;Yamamoto et al., 2017), this technique has rarely been applied to the evaluation of community assembly processes. In this study, the eDNA metabarcoding successfully reveals how the compositional variation changes in response to seasonal environmental shift, which illustrates its potential to be applied to community assembly studies. However, there remain several challenges for the technique to be widely applied. First, how accurately the eDNA signal reflects the actual species distribution is still an important issue to be explored. While its relatively short persistence and low diffusion rate promise its utility for detecting "snapshot" distribution Port et al., 2016;Thomsen et al., 2012;Yamamoto et al., 2016), whether the sampling is designed to meet the resolution must be examined carefully. For example, although our monthly sampling at sufficiently distant locations (>4 km apart) apparently enabled us to distinguish the temporal and spatial (horizontal) compositional differences, some of our sampling sites were very shallow, where the vertical mixing of eDNA potentially biases the detected distribution of species. Indeed, we found that the compositional difference between surface and bottom layers detected using eDNA metabarcoding analysis was smaller at shallower sites ( Figure S4) influenced not only by the presence of the species but also by the degradation and accumulation of the eDNA itself (Barnes & Turner, 2016). For example, the degradation of eDNA after released into water is faster at higher temperatures (Jo et al., 2019;Tsuji et al., 2017), and it has been assumed that eDNA accumulates in surface waters (Eichmiller et al., 2014;Moyer et al., 2014). We consider that these biases are not important here or at least do not change our main finding, because species richness in surface water was high in summer when the degradation is expected to be faster ( Figure S3), and the richness was generally lower in the surface layer where the eDNA accumulation is expected ( Figure S3).

| CON CLUS ION
In ecosystems, species assemblages vary according to environmental heterogeneity, but to a variable extent. We found that, in a coastal ecosystem, despite the larger environmental heterogeneity during the hypoxic season, species composition was more similar between bottom and surface layers, likely because the adverse effect of hypoxia was nonspecies specific. This result may seem to be contradictory to the current pervasive view that severe environment acts as strong species sorting. However, we speculate that these two contrasting conclusions are compatible: While the increasing environmental severity surely strengthens the filtering by sorting out vulnerable species, once it exceeds a certain level, the filtering would no longer work because too severe environment harms all the species in a similar way. Therefore, species sorting by environment may be greatest at the intermediate level of environmental severity, a hypothesis to be further confirmed in the future studies. Furthermore, by applying the eDNA metabarcoding analyses for examining community assembly processes for the first time to our best knowledge, this study highlights the potential applications of this promising technique across a wide range of disciplines of ecology.

ACK N OWLED G M ENTS
We thank Dr. Y. Miyazaki for helping with the filtering of fish species for the use of analyses with reference to records in Tokyo Bay. We thank two anonymous reviewers for their valuable comments. This work was supported by the Project for establishing a network of environmental and fisheries information.

CO N FLI C T O F I NTE R E S T
Authors declare no competing interests regarding this study.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this study are available in the Figshare repository (https://doi.org/10.6084/m9.figshare.19609347).