Habitat changes and population dynamics of fishes in a stream with forest floor degradation due to deer overconsumption in its catchment area

Forest degradation caused by deer overabundance has become a worldwide problem in recent decades. Overgrazing by deer not only affects terrestrial ecosystems but also spreads to aquatic ecosystems. Mass consumption of forest floor vegetation by deer creates denuded slopes and increases sediment inputs into adjacent rivers. In addition, rivers have upstream–downstream continuum structures, whereby the effects of degradation events in forests at upstream sites may spread to larger ecosystems downstream. However, few studies have examined the indirect effects of deer overabundance on downstream ecosystems. I examined the relationships between population dynamics of 13 fish species and habitat characteristics at a downstream site over the course of 11 years after forest floor degradation by deer overconsumption in a 36.5‐km2 catchment area of the Yura River in the Ashiu Research Forest, Japan, which is well‐protected from anthropogenic influences. During my 11 years of observation, characteristics of stream habitats changed from a predominantly coarse substrate to a fine substrate. I observed a remarkable decrease in one species (Tribolodon hakonensis) and increase in another species (Pseudogobio esocinus), and these changes were reasonably consistent with the increase or decrease in their preferred habitat types in the sampling site. This study showed long‐term habitat changes in a stream after forest floor degradation due to deer overconsumption in its catchment area and demonstrated that fish populations reacted to these changes. This study suggests that catchment‐level management, including forest ecosystem conservation, is necessary to solve fundamental problems in stream ecosystems.


| INTRODUCTION
In recent decades, forest degradation due to overabundant deer populations has become a problem worldwide (Côté, Rooney, Tremblay, Dussault, & Waller, 2004;Fuller & Gill, 2001). Deer overabundance not only affects terrestrial ecosystems but can also spread to aquatic ecosystems. Mass consumption of forest floor vegetation by deer creates denuded slopes and increases sediment inputs to nearby rivers (Sakai et al., 2012Biligetu et al., 2013). These increased sediment inputs affect the community structure of macroinvertebrate assemblages via changes in the physical characteristics of stream habitats (Sakai et al., 2012. In addition, rivers are continuous structures that transport waterborne materials downstream (Vannote, Minshall, Cummins, Sedell, & Cushing, 1980). It is therefore logical to hypothesize that when an overabundance of deer causes a large degree of forest degradation, the effects may spread through stream ecosystems to cause changes downstream. However, up to this point, very little research has focused on the indirect effects of deer overabundance on downstream ecosystems.
The effects of deer overabundance on downstream riverine ecosystems might be similar to those caused by increased sediment inputs at upstream sites due to human activities such as logging, agriculture, and urbanization. The loss of riparian vegetation due to human activity drastically increases the quantities of sediment and nutrients that enter streams from adjacent forests (Allan, Erickson, & Fay, 1997;Lowrance & Hendrickson, 1984;Osborne & Kovacic, 1993;Roth, Allan, & Erickson, 1996). Increased sediment yields upstream result in a lower proportion of coarse stone river bottom and a state shift from distinct habitats with deep pools and riffles to uniform channelized habitats (Allan et al., 1997;Roth et al., 1996). These changes usually decrease the quality of habitats for organisms that are adapted to local environmental characteristics such as the specific flow regime and habitat heterogeneities of the stream ecosystem in which they live (Fausch, Baxter, & Murakami, 2010;Gido, Dodds, & Eberle, 2010).
However, there are some difficulties inherent in determining the effect of deer overabundance on downstream ecosystems. If human activities and overgrazing by deer lead to similar downstream results, it will be difficult to separate the effects of humans and those of deer using an observational approach. In addition, the long-distance transportation of sediments from small upstream rivers to large rivers downstream is usually affected by seasonal and annual variations in discharge, and it may take place over the course of a year or longer (Church, 2002;Hamilton, 2012). Therefore, to examine the effects on downstream ecosystems of deer overabundance upstream, it is necessary first to identify a catchment area subject to deer overpopulation but well protected from human activity and then to perform long-term seasonal sampling and observations of fishes and habitat characteristics downstream.
In Japan, Sika deer (Cervus nippon) have become increasingly overabundant since the 1990s (Kaji, Miyaki, Saitoh, Ono, & Kaneko, 2000;Takatsuki, 2009). The distribution range of this species expanded by as much as 70% between 197970% between and 200270% between (Nakajima, 2007, and in some locations, the population increased almost 20-fold between 1988 and 2001 (Tsujino & Yumoto, 2004). The Ashiu Research Forest is located in the northern part of Kyoto Prefecture, central Japan, and contains the catchment area of the headwaters/upper reaches of the Yura River. This forest is managed by Kyoto University and is well protected from anthropogenic influences over a large area (https://data.ltereurope.net/deims/site/lter-eap-jp-58, accessed November 20, 2018). Although from the 1920s to the 1990s, the floor of the Ashiu forest was mostly covered with diverse herbs, small bamboos, shrubs, and woody seedlings (Yasuda & Nagamasu, 1995;Figure 1a), damage to forest floor vegetation caused by Sika deer browsing has become apparent since the late 1990s, and most of the forest floor had been denuded by the late 2000s (Fukuda & Takayanagi, 2008;Kato & Okuyama, 2004;Inoue et al., 2013;Tanaka, Takatsuki, & Takayanagi, 2008;Figure 1b,c).
In this study, I investigated the relationships between the population dynamics of 13 fish species and changes in habitat characteristics downstream from an area of forest floor degradation in the catchment area of the Yura River in the Ashiu Research Forest over a period of 11 years. I predicted that (a) increased sediment yields from upstream sites would decrease the proportion of coarse stones and increase that of fine stones in the streambed of a downstream site, as has been shown by previous studies (Allan et al., 1997;Gido et al., 2010;Roth et al., 1996); and (b) an increase or F I G U R E 1 Forest floor around headwaters (Kami-tani valley) in the Ashiu Research Forest in the summers of (a) 1998, (b) 2008, and (c) 2018 (left and center photos, by S. Shibata and T. Yoshioka, were provided by the Field Science Education and Research Center, Kyoto University) decrease in a certain habitat type would cause a corresponding increase or decrease in the population density of fish species that prefer that habitat. In this study, I also discuss the processes of the changes in habitat characteristics and fish populations that may have been induced by Sika deer, with reference to past climate and hydrological data. I set 22 line transects perpendicular to water flow along the channel at 10-m intervals at the sampling site, established six plots along each transect at regular intervals, and measured microhabitat environments (water depth, current velocity, and the proportions of sand, granule, pebble, cobble, boulder, and bedrock) at each plot. Data and detailed methods were archived in the database of the Japan Long Term Ecological Research Network (JaLTER; ERDP-2018-07; Nakagawa, 2019).

| Climate factors and hydrological events
The regional climate in the Ashiu Research Forest is warmtemperate with monsoonal effects. The annual precipitation is 2,257 mm, and the snow depth in winter averages~1 m at the sampling site and~2 m at the top of the catchment area. The principal flood events are attributable to snowmelt in spring, heavy precipitation in early summer, and typhoons in late summer.
Data on precipitation (daily total, mm) and snow depth (measured every morning, cm) at the sampling site (35 18 0 N, 135 43 0 E, 356 masl) from 1965 to 2018 were provided by the Ashiu Forest Research Station (https://data.lter-europe.net/ deims/site/lter-eap-jp-58, last accessed November 20, 2018). I obtained data on the water level of the Yura River at about 2 km downstream of the sampling site (135 41 0 E, 35 19 0 N, 279 masl) from the database of the Japan Ministry of Land, Infrastructure, Transport and Tourism for 1965-1988 and 2002-2018, and from that of Kyoto Prefecture for 1983-2010.

| Data preprocessing
I calculated the mean and SD of water depth and current velocity, and the mean proportion of each substrate type (sand, granule, pebble, cobble, boulder, and bedrock) on the streambed for each observation (Supporting Information).
I used the total number of individuals of the 13 most numerous fish species observed during each observation for analysis (Supporting Information). I counted the nocturnally active species Hemibarbus longrostris, Pseudogobio esocinus, Cobitis sp., Liobagrus reinii, Cottus pollux, and Odontobutis obscura during night observations and performed daytime counts of the other species. I did not use two fish species (Plecoglossus altivelis and Oncorhynchus masou) for statistical analysis despite their numerical dominance because they were artificially introduced by a local fishing society each year. Only individuals of 1 year or older were used for the analysis because of the high mortality rate and difficulty of field identification of fish larvae.

| Statistical analysis
I fitted temporal changes in fish populations and stream habitat characteristics to a state-space model (SSM; Commandeur & Koopman, 2007). Because measurements of habitat characteristics were not completed for all of the observations and the periods between observations were not entirely consistent, I estimated missing observation values when an observation did not have complete environmental data or when a period of non-observation was 2 months or longer (Supporting Information).
Process errors were assumed to follow a normal distribution in the SSM: where c is a constant to examine the directional change of habitat characteristics (Y) that was preliminarily assumed, μ t is the mean habitat state at time t, and σ μ is the SD of the mean habitat state. I assumed a normal error distribution for the mean and SD of water depth and current velocity as: where S t is the seasonal adjustment at time t and σ Y is the SD of observation errors. A six-step seasonal cycle was assumed as the seasonal adjustment in the SSM following Commandeur and Koopman (2007) (Supporting Information). I assumed a beta error distribution for the mean proportion of each substrate. Then, the process errors were assumed as: because the mean proportion of each substrate was described using the inverse-logit link function as: where φ Y is a divergence parameter. I described the process of temporal changes in fish populations using the Gompertz population model, a loglinear form of discrete logistic equation function (Fukaya et al., 2013;Royama, 1992;Turchin, 2003): where λ t is the mean state of a fish population at time t, σ λ is the SD of the mean state of a fish population, r s are the season-specific rate of increase of a fish population in season s at time t -1, and d are the density dependence, respectively (Fukaya et al., 2013;Royama, 1992;Turchin, 2003). I assumed season-specific parameters for the rate of increase of a fish population, because the focal fish species are known to have clear seasonality in breeding and growth (Kawanabe & Mizuno, 2001;Nakamura, 1969). As in the model of habitat changes, I estimated missing observation values for periods of 2 months or longer that held no observations and a six-step seasonal cycle was assumed (Supporting Information). In this model, the six-step seasonal cycle may represent a season-specific factor that is not related to population growth of fishes but effects on the count of fishes such as the finding rate of the individuals of fishes in a snorkeling observation in each season. I assumed a Poisson error distribution for the number of counted individuals of each fish species as: where Y t is the number of counted individuals of each fish species at time t, and α is the coefficient of the effect of habitat condition. I used the seasonal-adjusted level of the mean proportion of sand at time t (X t ), which was estimated from the SSM as the habitat condition variable, because I preliminarily assumed the increase of sand substrate in the sampling site after deer overabundance upstream, and this study aimed to examine the effects of the long-term trends in that change on fish population dynamics. In this model, the number of counted individuals of each fish species at time t (Y t ) was assumed to be determined by the past state of the fish population (Y t-1 ) and the present state of habitat characteristics estimated from the SSM (X t ). For missing observation values, I used a Gamma error distribution with a fixed shape parameter (=1) as: because an integer parameter has not been implemented in the statistical framework that I used (Stan Development Team, 2017). I conducted the analysis under a hierarchical Bayesian framework, using Markov Chain Monte Carlo methods via No-U-Turn sampling with the software STAN 2.17.0 (Stan Development Team, 2017), which is called from R 3.4.1 (R Development Core Team, 2017) using R package "rstan" (Stan Development Team, 2017), to obtain the posterior distributions of the parameters. I used the vague prior distribution Normal (0, 100) for the directional trend of habitat characteristics, the seasonal adjustment (S), the seasonspecific increasing rate (r) and density dependence of a given fish population (d), and the coefficient of the effect of habitat condition (α); and half-Cauchy (0, 5) for any variance and dispersion parameter terms in the examinations. I used four chains, each of which had 10,000 iterations including 5,000 burn-in with two thinning, resulting in 10,000 values for the posterior distribution of each parameter. Convergence of Markov chain Monte Carlo (MCMC) computing was checked using the "Rhat" and "n_eff" statistics (Rhat ≤ 1.1 and n_eff ≥ 100 for all parameters; Gelman et al., 2013). When the MCMC computing was not convergent, iterations and burn-in were increased up to 40,000 and 20,000, respectively, resulting in 40,000 values in the posterior distribution of each parameter. In addition, to confirm identifiability of the parameters with the vague prior distribution (S, r, d, and a), I calculated the overlap between prior and posterior distribution of the parameters (τ) as following Garrett and Zeger (2000), and checked all the parameters were τ ≤ 0.35 in each species (Garrett & Zeger, 2000;Gimenez et al., 2009).

| Changes in stream habitats
The mean proportion of sand clearly increased, while the mean proportion of boulder and the SD of current velocity '08 '09 '10 '11 '12 '13 '14 '15 '16 '17 '18 2007 2007 '08 '09 '10 '11 '12 '13 '14 '15 '16 '17  decreased during the observation period (Figure 2c,e,i). The sand proportion fluctuated around 0.05 from 2007 to early 2013, whereas the center of the fluctuation rose to 0.07 or higher after late 2014. Contrastingly, the mean proportion of boulder and the SD of current velocity fell from around 0.10 and 45 to 0.07 and 35, respectively, during the same period. The constants corresponding to the directional change in habitat characteristics estimated by the SSM analysis did not significantly deviate from zero (Table 1). In a comparison of microhabitats in the sampling site in June 2007 and June 2018, I observed a decrease in rapid current riffles and an expansion of stream beds covered by sand (Figure 3). The other habitat characteristics in the sampling site fluctuated but did not show any directional changes (Figure 2a,b,d,f-h,j).
During the sampling period, two of the three largest floods in the past 53 years occurred in November 2013 and October 2017, after large amounts of precipitation at the upper reaches of the Yura River (Supporting Information). Snow depth and spring floods relating to snow melt fluctuated annually, but without a clear trend (Supporting Information).

| Fish population dynamics
Identifiability of the parameters with the vague prior distribution were supported based on the τ estimate (Garrett & Zeger, 2000) for all the parameters in all the fish species with excepting Rhinchocypris oxycephalus, Pungtungia herzi, Hemibarbus longirostris, Niwaella delicata, and Cobitis sp. (Supporting Information). Among the eight examined fish species with the confirmation of the identifiability of the parameters, I observed a remarkable decrease in one species (Tribolodon hakonensis) and an increase in another (P. esocinus) between 2007 and 2018 ( Figure 4a,b), and the SSM indicated a significant negative effect of the change in habitat characteristics on the population dynamics of the first fish species and a significant positive effect on those of the second species (Credible Interval ≥ 95%; Table 2). From 2007 to 2011, T. hakonensis were constantly present in numbers of 100 or more, whereas after 2013, the annual maximum observation of this species was 50-100 individuals. By contrast, P. esocinus were rarely observed from 2007 to 2012, whereas after 2016 there were at least two individuals of this species in the majority of observations. Three fish species, L. reinii, Rhinogobius flumineus, and Nipponocypris temminckii, appeared to decrease, but this trend was unclear due to a large interannual fluctuations in populations, and the SSM did not indicate a significant effect of habitat changes on their population dynamics (Figure 4c-e; Table 2). Annual peaks of the number of individuals of these species demonstrated changes of twofold or more between neighboring years, but all of the largest and second-largest peaks in the populations of these fishes were observed before 2011. The other three fish species did not show any increasing or decreasing trends, and the SSM did not indicate a significant effect of habitat changes on their population dynamics (Figure 4f-h; Table 2).

| DISCUSSION
Temporal changes in the mean proportion of sand and boulder, as well as the SD of current velocity supported the prediction from the previous studies of the loss of riparian vegetation due to human activities (e.g., Allan et al., 1997;Roth et al., 1996). That is, the stream bed at the sampling site would shift from being predominantly coarse stones to predominantly fine stones, and the state of the stream would shift from distinct habitats with deep pools and riffles to uniform channelized habitats. Although I was unable to determine the exact timing of the habitat change due to the low frequency of observations during 2013-2016, the increase in sand substrate was not obvious before 2012. This pattern contrasts with the rapid response to forest degradation seen in headwater streams in the Ashiu Research Forest. Sakai et al. (2012Sakai et al. ( , 2013 reported that an increase in fine sediments in a streambed was clearly observable in first-order streams 2 or 3 years after the forest degradation event in 2006. This difference may be due to the time lag of downstream sediment transportation, which depends on the distance between upstream and downstream sites and seasonal and annual variations in water discharge, including flood disturbance events (Church, 2002;Hamilton, 2012). In small tributaries (i.e., first-or second-order streams) of mountain streams, sediments are directly supplied from neighboring slopes and are transported downstream; this is the main process by which higher order streams (i.e., ≥third-order streams) acquire fine sediments (Church, 2002). In addition, sediment transport to downstream sites mainly occurs during large flood events. Large proportions of these fine sediments arrive and are stored in riparian sites of river channels, thereafter being gradually supplied to streambeds by ordinal rains and floods over a period of months or years. The nonsignificance of the directional trend in the SSM analysis despite the obvious change in habitat characteristics at the sampling site may be explained by many missing observation values and the dependence of sediment transportation on flooding. I could not conduct frequent samplings from 2013 to 2016, and that timing seems consist with that of the increase of the proportion of sand. In addition, large floods occurred in September 2013 and October 2017 in the sampling site. I suggest that a large amount of fine sediment was transported from upstream sites to the riparian area of the sampling site at those times and that the stored sediments might have supplied fine substrate materials to the streambed by the process explained above. However, that directionalchanging event of habitat might not reflected to the directional trend parameter (c) in the SSM due to the large variance of missing observation values.
The patterns of habitat changes and population dynamics of fish were reasonably consistent for at least two species. T. hakonensis use pebble-cobble-dominated streambeds with the effect of subsurface water for spawning (Kawanabe & Mizuno, 2001;Nakamura, 1969) and prefer pool habitats with shelters such as accumulations of large rocks or deep crevices in bedrock (Katano, Aonuma, & Matsubara, 2001). The decrease in this species may have been caused by the decrease in its preferred habitats at the sampling site. For example, shelter or spawning habitats preferred by this species might have been buried by sand supplied from the upstream catchment area. In contrast, P. esocinus increased between 2007 and 2018. This benthic species is a vacuum feeder with a characteristic mouth that can elongate in the anteroventral direction. It forages macroinvertebrates in fine substrates, and often buries itself in sand when threatened (Kawanabe & Mizuno, 2001;Nakamura, 1969). In the Yura River in 2008, P. esocinus mainly inhabited sand-dominated sites downstream of the sampling site (Nakagawa, 2014). The increase in fine sediment substrates at the sampling site may have expanded the range of distribution of P. esocinus from downstream sites. These consistent patterns support the conclusion that deer overconsumption near the catchment area of a river can have indirect effects on the population dynamics of fish downstream.
Several fish species (e.g., C. pollux and Cobitis sp.) did not clearly decrease or increase at the sampling site during the 11 years of observation despite their preference for or aversion to coarse substrate habitats (Kawanabe & Mizuno, 2001). The population dynamics of these species may have been driven by factors other than those examined in this study, for example, direct mortality from flood disturbances (Jensen & Johnsen, 1999;Kamp, 2017) and interspecific competition (Harwood, Mecalfe, Armstrong, & Griffiths, 2001;Zwol, Neff, & Wilson, 2012). Additionally, the response of a fish population to habitat changes is usually nonlinear and has a species-specific threshold (e.g., Gido et al., 2010;Perkin & Gido, 2011). In this study, the population density of T. hakonensis clearly decreased from 2009 to 2012, remaining constant at a low density thereafter, 2007 '08 '09 '10 '11 '12 '13 '14 '15 '16 '17 '18 2007 '08 '09 '10 '11 '12 '13 '14 '15 '16 '17 '18 2007 '08 '09 '10 '11 '12 '13 '14 '15 '16 '17 '18 2007 '08 '09 '10 '11 '12 '13 '14 '15 '16 '17 '18 Number  whereas P. esocinus noticeably increased in the most recent few years. Fish species that did not show any obvious population changes during the period of this study may merit further monitoring over a longer period of time. This study showed long-term habitat changes in a stream with a large catchment area protected from anthropogenic influence and demonstrated that fish populations significantly reacted to these habitat changes after forest floor degradation due to deer overconsumption in the catchment area of the stream. Ecosystem linkage greatly affects population productivity and consumer-resource dynamics in the recipient ecosystem of a resource subsidy (Polis, Anderson, & Holt, 1997), and habitat connectivity has important roles in the stability of local populations and the dynamics and structure of local communities (Leibold et al., 2004). The positive roles of these ecosystem/habitat linkages have been widely discussed (e.g., Baxter, Fausch, & Saunders, 2005;Nakano & Murakami, 2001;Terui et al., 2018b). However, these linkages also propagate unfavorable effects, such as the dispersal of alien species (Atobe, Osada, Takeda, Kuroe, & Miyashita, 2014;Rahel, 2007), pathogens (Jules, Kauffman, Ritts, & Carroll, 2002), and pollutants (Hashimoto et al., 2013;Lebowitz, 2003) and these indirect effects (Paetzold, Smith, Warren, & Maltby, 2011;Terui et al. 2018a). It is well known that the top-down effects of large herbivores greatly change terrestrial ecosystems (Côté et al., 2004;Fuller & Gill, 2001;Paine, 2000). This study presents a novel finding: A large degree of environmental change due to top-down effects of large terrestrial herbivores in one ecosystem can harm adjacent ones via ecosystem linkages.
Additionally, this study suggests that catchment-level management, including conservation of forest ecosystems, is necessary to solve fundamental problems in stream conservation. At present, the concerns of governments, forest managers and scientists regarding Sika deer have been limited to commercial damage and/or biodiversity loss on land (e.g., Kaji et al., 2000;Takatsuki, 2009). However, the problems caused by Sika deer might not be restricted to terrestrial issues in situations wherein stream habitats may be affected via ecosystem linkages, as seen in this study. In that case, Sika deer cannot simply be dismissed by fisheries managers and stream scientists as "someone else's problem." The types of habitat changes recorded in the observations in this study are likely to affect organisms other than fishes, such as algae and aquatic insects, and even to have an effect on the total structures of stream food webs and ecosystem processes. Although researchers, a local government, and local peoples have attempted to control the deer population in Ashiu Research Forest, there have been no signs of recovery in the forest floor vegetation (Inoue et al., 2013, Figure 1c). It seems likely that the changes in stream habitats seen in this study may progress for at least several more years. In addition, degradation of forest floor vegetation has expanded to large areas around the research forest (Takatsuki, 2009). Deer overabundance may therefore be affecting larger aquatic ecosystems than the sampling site. In this study, I only observed changes in fishes and their habitats after forest floor degradation. In future works, I hope to reexamine the effects of deer based on long-term monitoring data, including fish population dynamics and habitat changes after recovery of forest floor vegetation in the river catchment area.
ACKNOWLEDGMENTS I thank S. Shibata, T. Yoshioka, and the Field Science Education and Research Center of Kyoto University for providing the old photos of the Ashiu forest. I thank the staffs of Ashiu Forest Research Station, the staff of Miyama Gyokyo (Fisheries Cooperative of Miyama, Yura River), and local people in the Ashiu for kind cooperation to my research. This study was supported by JSPS KAKENHI Grant Number JP19K15857 and Global COE Program A06 "Formation of a Strategic Base for Biodiversity and Evolutionary Research: from Genome to Ecosystem" from the Ministry of Education, Culture, Sports and Technology, Japan.

AUTHOR CONTRIBUTIONS
H.N. contributed all the parts of this study, planning, samplings, analysis, and writing a paper.

ETHICS STATEMENT
Material fishes never have been directly caught in this study.