Factors explaining the yearly changes in minimum bottom dissolved oxygen concentrations in Lake Biwa, a warm monomictic lake

Vertical profiles of dissolved oxygen (DO) and water temperature (WT) measured bi-monthly for 36 years (1980–2015) near the deepest part of a warm monomictic lake were analyzed with special reference to yearly minimum DO at bottom (DOmin). DOmin changed yearly (3.0 ± 1.2 mg l−1) and significant differences in DOmin were not observed between Period I (1980–1993; cooler and worse in water quality) and Period II (1994–2015; warmer and better in water quality). This unclear trend in DOmin was probably due to the offsetting influences between warming induced by global warming and oligotrophication attempted by local governments etc. for the study period. DOmin was positively correlated with disturbance time (timing of last cold water intrusion observed from Mar to Aug), which could be related to the start of DO depletion at bottom. Thus, the linear model using this parameter could predict yearly DOmin fairly well for the entire study period (r2 = 0.60). In addition, DOmin and time of disturbance were correlated negatively with water density at bottom in Jan and positively with water density equilibrated to air temperature (AT) in Mar. Higher lake water density after full depth mixing advances the disturbance time. In contrast, lower AT in Mar and/or higher density of influent water after Mar delays the time likely due to the larger amount of snowfall in the watershed. Further, DOmin was positively correlated with maximum wind velocity in Sep which probably induced the recovery of DO. Multiple-regression models to predict DOmin using these meteorological and water quality parameters were developed (r2 ≥ 0.38, worse performances than the model using disturbance time) to forecast future trends of DOmin through global warming and/or climate change. Significant influences of water or sediment oxygen demands on DOmin were not detected. We also discuss the applicability of the proposed models.

The dynamics of dissolved oxygen (DO) distribution in inland waters are fundamental to the development of an understanding of the distribution, behavior, and growth of aquatic organisms 1 . The DO state reflects the health of the lake environment. Recently, there have been concerns about the effects of global warming on the state of lake DOs, particularly bottom DO, through physical, chemical and biological processes 2 . As a result, the Japanese government has started to use bottom DO conditions as an environmental standard in lake and coastal water bodies 3 .
In European and North American lakes, particularly warm monomictic or dimictic lakes, depletion of DO in the hypolimnia during stratification has been observed and analyzed for more than 100 years 4,5 . The degree of such hypolimnetic DO depletion can be affected by the duration of the stratification period and the DO depletion rate 5 . In general, the DO deficit is expected to grow owing to eutrophication and warming 1 . Similarly, hypolimnetic waters have been reported to become anoxic in deep tropical lakes 6 . In response to these changes, a number of diverse approaches have been adopted to investigate the temporal dynamics and spatial patterns of hypoxia in European lakes and coastal regions 2 .
In Lake Biwa (Fig. 1), the largest lake and most important freshwater resource in Japan, much scientific and social attention has been paid to bottom DO problems in recent decades. The north basin of Lake Biwa is itself a SCiEntifiC REpoRTs | (2019) 9:298 | DOI: 10.1038/s41598-018-36533-7 warm monomictic lake. Using limnological observation data up to 1971, Naka 7 showed that the amount of DO in deep water before the season of full-depth mixing has been decreasing due to a progressing eutrophication. Sohrin et al. 8 analyzed long-term changes in bottom water, and showed that, since 1999, yearly minimum DO concentrations <50 μmol kg −1 have started to occur frequently at around a depth of 90 m. Although they considered that these changes were most likely the result of global warming (S- Fig. 1) and local eutrophication, they did not specifically analyze the cause(s). Fushimi 9 suggested that a trend of decreasing snowfall in the catchment might be contributing to the decrease of DO in the deep layer. Based on an analysis of the vertical profiles of DO and influential factors, Tsujimura et al. 10 and Jiao et al. 11 reported that the DO conditions in the hypolimnion during the stratification period were  affected by the occurrence of violent lake-water mixing (possibly due to meteorological events, e.g., typhoons), the species and amount of phytoplankton in the epilimnion in summer, and water temperature and DO concentration in spring. They also investigated the impact of a low DO condition on benthic organisms. Okamoto 12 described the changes in bottom DO concentrations near the deepest part of the lake during specific years in great detail, and discussed the influencing factors. These analyses all considered that the yearly minimum DO might be related to the beginning of stratification, time of overturn, and rate of decrease of DO (DO decrease rate) under the stratification conditions. However, no report has quantitatively explained and/or predicted the yearly minimum DO in bottom water in recent decades.
The purpose of the present study was therefore to elucidate the factors affecting yearly minimum DO concentrations in bottom water using bi-monthly data from the last 36 years on the vertical distribution of water temperature and DO at the center of Lake Biwa. Much attention was paid to meteorological factors, since it is widely anticipated that such factors will reveal the influences of global warming and/or climate variation on lake DO conditions. Such information would be useful not only for building a detailed simulation model describing the water cycle, and physical, chemical and biological processes in the watershed and lake, but also when considering countermeasures for low DO problems.

Results
Seasonal and yearly change in DO and WT. The observed changes in DO and water temperature (WT) at the lake bottom are shown in S- Fig. 2. Seasonal cycles in DO are clearly indicated, and gradually increasing oscillations of several-year periods overlapping on seasonal cycles were observed in WT. Similar oscillations in WT were also observed at St. 2 8 . Examples of seasonal changes in vertical profiles of DO and WT for specific years are shown in Fig. 2 and S- Fig. 3. Generally, DO at bottom decreased monotonously from the start to the end of the DO depletion period, but occasional small disturbances (recoveries of DO) were also observed, e.g., Jul 2008 in Fig. 2 (2). This disturbance was possibly caused by a horizontal movement of bottom waters, but the details are unknown.
As can be seen in Fig. 2 and S- Fig. 3, DO at bottom decreased substantially in 1987, 2002 and 2008, while the decreases in 1983, 1990 and 2014 were smaller. The starting times of DO decrease occurred earlier in the former years compared to the latter. WT changes appeared roughly similar in all six of these years, but the minimum WT at bottom (indicated by arrows in the figures) was observed sooner in the former years than in the latter. This timing of minimum WT at bottom seemed to influence the DO decrease at bottom. Thus, we call this the disturbance time, and analyze it below. Changes in the yearly minimum DO at bottom (3.0 ± 1.2 mg l −1 : Fig. 3), the DO minimum date (S- Fig. 4), the WT of the disturbed water, and the date of disturbance did not show clear trends or cyclic oscillations (Table 1).

Relations between the yearly minimum DO and influencing factors.
Results of the correlation analysis between yearly minimum DO concentrations and yearly values of the factors possibly influencing it were as follows (*P < 0.05, **P < 0.01). At first, the yearly minimum DO (DOmin) was not significantly correlated with specific WT or DO values, e.g., minimum WT at bottom, WT at bottom when the minimum WT at the surface was observed, DO when the minimum WT at the surface was observed, and WT when the minimum DO was observed. DOmin showed non-significant correlations with the dates of the unusual WT or DO values, e.g., the dates when minimum WT at the surface, minimum WT at bottom, and DOmin were observed. In addition, non-significant correlations were found between DOmin and the periods indicating DO depletion duration, e.g., the periods from the time of minimum WT at the surface or bottom to the time when DOmin was detected (the time of DOmin). With regard to meteorological factors in the spring, each of the following showed a higher correlation with DOmin in a spring month compared to other months: air temperature (AT) in Mar, precipitation (PR) in Apr, and sunshine hours (SH) in Apr at Imazu. However, none of these correlations were statistically significant. The snowfall amount (SA) at Imazu during winter showed a positive correlation with DOmin, but this relation was also not significant. As to the values in fall, the correlation coefficients with AT in Sep and monthly maximum wind velocity averaged for 10 minutes (max WV) at Hikone in Sep (r = 0.34*; r: correlation coefficient) were high. The number of cold days did not have a significant influence on DOmin. Regarding WT at bottom, this parameter was significantly correlated with DOmin at the Jan 1 st survey (r = 0.43**) and Jan 2 nd survey (r = 0.33*), but the correlations with WT in the fall, e.g., WT at the surface at the Oct 1 st survey were non-significant. In addition, a significant correlation was found between DOmin and the DO decrease rate from Apr 1 st to Sep 2 nd survey (r = −0.49**).
Because the density of water (ρ) (a function of water temperature as shown in the Methods) has a more direct influence on the mixing process in lakes than temperature, we analyzed the correlation between DOmin and water density. Compared to temperature, higher correlations were found between DOmin and ρ as a function of WT at the Jan 1 st survey (r = −0.43**), and DOmin and ρ of water equilibrated to AT in Mar. As suggested in the  The correlation between DOmin and yearly averaged water quality-e.g., biochemical oxygen demand (BOD), chlorophyll a, total phosphorus (TP)-was not significant. The rate of decrease in DO from the Apr 1 st to Sep 2 nd survey was significantly correlated with the rates of decrease for other long periods (i.e., periods of more than 3 months), but not with the rates of decrease for short periods (less than 2 months). The DO decrease rate from the Apr 1 st to Sep 2 nd survey showed positive correlations with ATs, but those correlations were not significant.
Water quality and air temperature trends during the survey period. In the north basin of Lake Biwa, an oligotrophication trend was observed for the period 1980-2015, as shown in S- Fig. 6 13 . In addition, the shift analysis indicated clear changes between 1993 and 1994 for BOD and TP and between 1999 and 2000 for chlorophyll a. Thus, we divided the whole period into two parts: Period I from 1980-1993 and Period II from 1994-2015. There were significant differences between the two periods in the concentrations of BOD, chlorophyll a, and TP, averaged surface WT from Aug-Oct., yearly minimum bottom WT, average of bottom WT from Aug-Sep, WT of the disturbed water, average of AT from Mar to Oct, and number of cold days from Oct 1 to Nov 15 (Table 1). In Period II, the concentrations of all the parameters related to water quality (BOD, chlorophyll a, TP) were lower except the nitrogen level. WTs except WT at DOmin and ATs were higher in this period. In general, lake water quality was better and the lake environment was warmer in Period II than in Period I.
Annual averaged air temperature increased at a rate of 0.039 °C y −1 (the regression coefficient in the linear model with year) during the period from 1980-2015. The increase rates were different month-to-month (S- Fig. 2 (2)). Significantly larger increases were observed from May-Oct, while smaller, non-significant increases were observed in other months.

Prediction of DO minimum in the bottom layer.
First, models were constructed using the variable (i.e., the disturbance time) that had the highest correlation coefficient with DOmin in this study. Its linear, quadratic and cubic equations were compared based on the adjusted r 2 values and then the following linear one was selected (n = 36, r 2 = 0.60**, adjusted r 2 = 0.58, RMSE = 0.74 mg l −1 ): × − . DOmin 0 027 Tdisturbance 1 484 (1) where Tdisturbance is disturbance time (Julian date). The characteristics of yearly changes were well described, but the amplitudes of the changes were somewhat smaller than the observed values due to the theoretical characteristics of regression model (Fig. 3). Next, the multiple regression models to predict DOmin were examined using meteorological and water quality parameters. Based on the correlation analysis shown above, we compared the models with combinations of the following parameters as independent variables: WTs equilibrated with AT at Imazu in Feb and Mar, water densities of these WTs, bottom WTs during the first and second periods of Jan, bottom WTs during the first and second periods of Feb, water densities of these WTs, total sunshine hours at Imazu in Mar and Apr, precipitation at Imazu in Apr, monthly max WV averaged for 10 minutes in Sep at Hikone, and ATs at Imazu in Jul, Aug, Sep and Oct. Then, the following model was chosen to give the highest adjusted r 2 for all data (1980-2015) (n = 36, r 2 = 0.38**, adjusted r 2 = 0.30, RMSE = 0.92 mg l −1 : Figs 3 and 6): where DOmin is mg l −1 , ρ′: (ρ − 1) × 10 3 , ρ is the water density corresponding to the water temperature (g cm −3 ), WTJ1 is the bottom WT during the first period of Jan, WT(AT3) is the WT equilibrated with AT at Imazu in Mar, SH4 is the total sunshine hours at Imazu in Apr (hr), and W9 is the monthly max WV averaged for 10 minutes in Sep at Hikone (m s −1 ). Similar to eq. (1), this model gave somewhat smaller variations compared with the observed ones due to the rather small regression coefficient, and the discrepancies were larger in the latter period compared with the former. To clarify the differences between the periods, we applied multiple regression analysis separately for Periods I and II. The models obtained for the highest adjusted r 2 were eq. In eq. (3), the same parameters were chosen as in eq. (2). In eq. (4), a new parameter AT7 (air temperature in Jul) (°C) was used instead of SH4. The correlation between the predicted and observed DOmin (S- Fig. 9; combined eqs (3) and (4); r 2 = 0.44**, adjusted r 2 = 0.30, RMSE = 0.87 mg l −1 ) was better than the correlation shown in Fig. 6, and higher amplitudes of change were obtained. Particularly in Period I, better agreement was achieved.

Discussion
Characteristics of long-term trends in the DO minimum in bottom water. Sohrin et al. 8 reported that the bottom DO at St. 2 has been decreasing since 1959. They considered that the increases in bottom WT and nutrient concentrations (NO 3 -N, PO 4 -P, and Si(OH) 4 ) were responsible for the decrease in DO. The relationships between yearly minimum DO concentrations at St. 1 and St. 2 appeared to be rather scattered, but were significantly correlated (S- Fig. 7).
In this study, we could not identify a clear trend in the bottom DO concentration ( Fig. 3; Table 1). The difference in the periods investigated likely affected the difference in bottom DO trends and main factors influencing the DOmin. We therefore statistically analyzed the data from Sohrin et al. 8 14 . Due to such attempts, many parameters of water quality began to improve after the late 1990s 13,14 . There were gradual increases in WT and NO 3 -N concentrations after the late 1990s, but the changes in bottom DO after 1980 were not significant (S- Table 1). In addition, there were significant correlations between DOmin and minimum bottom WT (r = −0.31*), and between DOmin and yearly averaged NO 3 -N (r = −0.55**) for data between 1963 to 2012, but the correlations between DOmin and minimum bottom WT and between DOmin and yearly averaged NO 3 -N were not significant for the period between 1980 and 2012. In short, the study period (1980-2015) exhibited overall trends of gradual oligotrophication and warming, as discussed in the previous section and shown in Table 1. The unclear trend in yearly DOmin is probably attributable to the offset between suppression of DO deficit by oligotrophication and its intensification by warming during our study period.

Mechanism underlying the DO minimum in bottom water. Theoretically, DOmin concentrations
at the bottom of warm monomictic and dimictic lakes are determined by the duration of the stratified period, i.e., from its start to its end, and the DO depletion rate (e.g., Walker 5 ). However, significant correlations between DOmin and DO depletion periods were not obtained in this study. The rates of decrease of DO over long periods (more than three months) were significantly negatively correlated with DOmin. However, this correlation was merely due to different characteristics of the same phenomena, because these DO decrease rates were calculated as the amount of decrease in DO (closely related to DOmin in the case of long period) divided by the duration of the period.
In contrast, the decrease rates for shorter periods (less than two months) did not show a clear influence on DOmin, indicating the scattering of the DO decrease rates during the stratification period. Tsujimura et al. 10 reported a significant correlation between the average volumes of phytoplankton cells from Mar to Sep and DOmin (r = −0.57*). In this study, we did not include such a parameter explicitly due to the imperfect availability of data over the whole study period. In eq. (4), AT in Jul was selected as one of the independent variables for the prediction of DOmin in Period II, suggesting that this parameter is related to the DO decrease rate.
Further, the sediment oxygen demand (SOD) might contribute to the determination of the value of DOmin 4 . There has been no report on SOD in the north basin of Lake Biwa. Gelda et al. 15 built a sediment diagenesis model expressing carbon and oxygen dynamics in water and sediment and then predicted the change in SOD over a 30-year period in Onondaga Lake. However, yearly estimation of SOD is quite difficult without deliberate observation or a dynamic model. Rucinski et al. 16 expressed SOD as a function of TP and then predicted the hypoxia in Lake Erie. Adopting their approach, we added TP to eqs (1)-(4) as an independent variable, but this yield no improvement in the adjusted r 2 . In addition, regression equations similar to eq. (1) were obtained separately for Periods I and II. There was no significant difference between these equations, indicating that the difference in SOD between the periods was not explicitly expressed in these equations. This may have been because there was an insignificant difference in SOD between the two periods, or because SOD did not have a significant influence on DO change during the isolation period of bottom water after the last intrusion of cold water. In any case, in order to accurately forecast the future hypolimnetic DO conditions in this lake, it would be necessary to track and predict the long-term change in oxygen consumption rates in water and sediment.
In Lake Biwa, the behavior of the hypolimnetic DO has been discussed from the perspective of the timing of complete oxygenation of the bottom water, the WT at that time, SA in the basin, events such as strong winds (e.g., typhoon), cold days in autumn, and phytoplankton species and their amounts in summer [9][10][11][12]17 . In this study, we analyzed the relations between DOmin and these influencing factors, and found that a limited number of factors (WT at the Jan 1 st survey, WT at the Jan 2 nd survey, max WV in Sep, water density as a function of AT in Mar), were significantly influential, but not predominant expressed by the relatively low regression coefficients in eqs (2)-(4).
However, the quite high correlation (Fig. 4: r = 0.77**) with the disturbance time, which determines the start of DO depletion at bottom, seemed likely to have affected the DOmin at bottom. In this sense, the timing of the last intrusion of cold water, when DO is plentiful (DO of the disturbed water: 9.0 ± 0.8 mg l −1 ), would be the start of DO depletion at bottom. The rather small variations in DO of the disturbed water indicated that the DO state was relatively constant at the start. Yoshiyama et al. 17 suggested the importance of gravity currents, as well as the contribution of convective mixing to the delivery of oxygen to the profundal zone of the north basin of Lake Biwa. Because a high correlation was observed between AT in Mar and SA at Imazu (r = −0.61**), low AT in Mar (i.e., high water density equilibrated to AT in Mar) could delay the disturbance time by melting the snow in the watershed more slowly and increasing the density of the influent water after Mar. It would be reasonable to consider that both the water density state in lakes after full-depth mixing and the water density of influent water after Mar affects the disturbance time, and therefore, the density difference between them could be a good predictor of disturbance time (S- Fig. 5 and Fig. 5) and DOmin (r = 0.49**). This density difference was a better predictor of the DOmin than the difference of AT (r = 0.39*). Importantly, higher lake water density after full-depth mixing advances this disturbance time, and vice versa. In contrast, higher density of influent water after Mar due to lower AT in Mar delays the time and vice versa.
Surprisingly, there was a positive, albeit non-significant, correlation between DOmin and its timing (i.e., the end of DO depletion). This result appears to indicate that, with respect to DOmin, the starting time of DO depletion is more important than the time of ending time. Because dissolved and particulate organic matter decompose logarithmically in lakes 18,19 , the initial processes (e.g., supply of organic matter to the bottom water) were suggested to have a great influence on the decrease compared with the subsequent processes. Only max WV at Hikone in Sep significantly raised the DOmin. Jiao et al. 11 showed that the wind impact index (defined as the square of wind velocity divided by Schmidt's Stability Index) affected the recovery of DO depletion. Okamoto 12 found that a typhoon passing near the lake accompanied by strong winds in autumn affected DO conditions at the bottom. In addition, wind velocity at Imazu had a positive, although non-significant, correlation with DOmin, suggesting that the wind in flat areas more strongly affects the mixing in the lake compared with the winds in mountainous areas. Meanwhile, there was a close relationship between AT for Sep to Dec in the preceding year and WT in the Jan 1 st survey (S- Fig. 8). A seasonally nonuniform increase in AT during our investigated period (e.g., a higher increase from Sep to Dec compared with Mar) could have resulted in a negligible delay in the start time of DO depletion in spite of the significant rise of annual mean AT. Thus, the seasonality of AT increase could affect the change in DOmin. For example, a change in WTJ1 from 7.7 °C (average for 1980-2015) to 8.7 °C would bring about an increase of 0.45 mg l −1 DOmin using eq. (2), assuming that other meteorological conditions were not varied. Future studies considering such phenomena will be necessary to understand and predict the influence of global warming and/or climate change on lake environments. Differences in the model predictability between Period I and II. The model agreement in Period I was better than that in Period II. Similar to eq. (1), linear models predicting DOmin were determined separately for Periods I and II. The RMSE (0.62 mg l −1 ) of the obtained model for Period I (Y = 0.0225 × −0.949 where X is the disturbance time and Y is the DOmin) was smaller than that (0.76 mg l −1 ) for Period II (Y = 0.0306 × −1.929), indicating a closer relationship with the disturbance time in Period I compared with Period II.
In the south basin of Lake Biwa, submerged plants hardly existed from 1960 to 1993, but they began to recover after the serious shortage of water in the summer of 1994 20 . Along with this recovery, there was an improvement of the transparency. It would be possible, albeit unlikely, that this regime shift occurring in the south basin led directly to the shift of water quality observed in the north basin, because there were large differences in lake morphology and water quality between the basins.
Kishimoto et al. 21 reported that the average cell size of the phytoplankton community in the north basin decreased from 269 μm 3 cell −1 in the 1980s to 56 μm 3 cell −1 in the 2000s; in particular, the proportions of small nanoplankton (mainly Cyanophyceae) increased in summer. After the year 2000, the yearly variations in the summer phytoplankton composition increased. Yamada et al. 22 indicated that the biodegradability of dissolved organic matter released from phytoplankton differed species by species based on experiments using Microcystis aeruginosa, Staurastrum dorsidentiferum, and Cryptomonas ovata under the water conditions existing in the north basin. These facts indicate that the larger yearly variability in the phytoplankton community for Period II compared to Period I resulted in less agreement in the DOmin prediction models. Therefore, biological information is necessary to improve the models, particularly in Period II.
Applicability of the DOmin-predicting model and future study. The prediction model in eq. (1) has a fairly simple form and shows a good performance; thus a bi-monthly monitoring of bottom WT and DO can provide an early warning signal of low DOmin several months before the event. However, this model cannot forecast the future trends of DOmin through global warming and/or climate change. In contrast, the prediction models in eqs (2)-(4) are less accurate, but they can roughly forecast the future trends in DOmin corresponding to climate change. These models could be used to evaluate more elaborate simulation models describing the water cycle and the physical, chemical and biological processes in watersheds and lakes. However, we should keep in mind that eqs (3) and (4) are essentially disadvantageous for the prediction of future lake conditions. At any rate, additional efforts will be needed to enhance the accuracy of these prediction models. Such efforts would be expected to include the collection of DO measurements with higher temporal and spatial resolution, meteorological data at more sites, and information on DO demand by water and sediments.

Methods
Lake Biwa. Lake Biwa, located at the center of Honshu Island, is the largest lake (area 670 km 2 ; maximum depth, 104 m) in Japan (Fig. 1). One of several ancient lakes 23 , it contains more than 50 endemic species (benthic fauna and deep-living fish), any of which may be vulnerable to low oxygen concentrations in profundal habitats 17 . Lake Biwa has two basins-one to the south and one to the north. The south basin is shallow, with a mean depth of 4.0 m, while the north basin is deep, with a mean depth of 43 m. Judging from the trophic boundaries proposed by Carlson and Simpson 24 , the current trophic level of the north basin is between oligotrophic and mesotrophic in terms of total phosphorus, and mesotrophic in terms of chlorophyll a. This basin is warm monomictic; a clear thermocline is formed in summer and fall, and full depth mixings occur in late fall and spring.

Survey and database.
A bi-monthly limnological survey (1 st and 2 nd surveys in the respective months) has been conducted at the center of Imazu-oki (St. 1 in Fig. 1: 35°23′41″N, 136°07′57″E, ~90 m depth) from Apr 1979 to the present by LBERI (Lake Biwa Environmental Research Institute). WT and DO have been measured at 0.5, 5, 10, 15, 20, 30, 40, 60, 80 m below the surface and at 1 m from the bottom sediment (hereafter, at bottom) using a Hydrolab Quanta Multi-Probe Meter (Hydrolab, Loveland, CO) calibrated by JIS k02102-32. During the period from 1980-2001, WT and DO 1-m above the sediment were measured after the system was laid on the sediment and then raised 1 m. After 2002, a system with a depth sensor was used, and values 1-m above the bottom were determined after data processing.
The analysis in this study used the data from 1980-2015, thus, 36 years of data were analyzed. During the periods when the observations were more frequent than bi-monthly, we excluded any irregularly timed observations (there were only a few such occasions over the 36-year observation period). There were also a few periods when the observations were less than bi-monthly, and approximately 10 occasions when no measurements were conducted at several depths. During these periods, we did not interpolate the values except for special analyses (i.e., determination of the DO decrease rates, as shown below).
As for other water quality values-e.g., BOD, chlorophyll a, TP, total nitrogen (TN), and nitrate nitrogen (NO 3 -N)-the values averaged for the Japanese fiscal year (Apr to Mar) and in the north basin were used 13 Fig. 1 measured by the Shiga Prefectural Fisheries Experiment Station were also used for comparison 8 . Daily data on meteorological factors, including AT, PR, SH, SA, max WV at Imazu, and max WV at Hikone (Fig. 1) were used for the analysis 25 . The changes in monthly averaged AT at Imazu are shown in S- Fig. 1. Water density. Based on the database 26 , the following quadratic function was used to calculate the water density (y: g cm −3 ) from water temperature (x: 0 to 10 °C). = − . + . + . y 00000077912 x 00000630498 x 0 9998445290 2 Indices for characterizing WT and DO changes and meteorological conditions. The following indices for water quality were determined and used for further analysis: (1) the yearly minimum DO at bottom (sometimes including the data from Jan in the following year depending on DO seasonal change) and its date (all dates are Julian); (2) the DO decrease rate (several periods); (3) the WT and date (referred to as disturbance time) of the bottom water disturbance event (after Apr; WT at bottom layer was lower than before and after the event; more than 8 mg l −1 , i.e., nearly DO-saturated water); (4) the yearly minimum WT at the surface and bottom waters and their dates. As a meteorological index, we calculated the number of days for several different periods during which the daily minimum AT was colder than the critical WT (number of cold days). The critical WT was that at the bottom when the yearly DO minimum was observed.
Statistical methods. Statistical analyses were used to determine the correlation, the differences between the means (t-test) and the multiple regression for significance at the level of 0.05 (*) or 0.01 (**) (Excel-Statistics for 2016: BellCurve Social Survey Research Information Co., Tokyo). The values of adjusted r 2 were used to select the independent parameters in multiple regression analysis. The root mean square errors (RMSE) were also calculated to evaluate the obtained prediction models. Here, r 2 is the square of the correlation coefficient (determination coefficient), N is the number of data, p is the number of independent variables, and X pred,i and X meas,i are the values of predicted and measured values, respectively.
The sequential t-test regime shift detection software ver. 3.2 27 was used to determine shifts in the time series (significance level = 0.1; cut-off length = 10; Huber's weight parameter = 1).