Linking Stoichiometric Organic Carbon–Nitrogen Relationships to planktonic Cyanobacteria and Subsurface Methane Maximum in Deep Freshwater Lakes

Our understanding of the source of methane (CH4) in freshwater ecosystems is being revised because CH4 production in oxic water columns, a hitherto inconceivable process of methanogenesis, has been discovered for lake ecosystems. The present study surveyed nine Japanese deep freshwater lakes to show the pattern and mechanisms of such aerobic CH4 production and subsurface methane maximum (SMM) formation. The field survey observed the development of SMM around the metalimnion in all the study lakes. Generalized linear model (GLM) analyses showed a strong negative nonlinear relationship between dissolved organic carbon (DOC) and dissolved inorganic nitrogen (DIN), as well as a similar curvilinear relationship between DIN and dissolved CH4, suggesting that the availability of organic carbon controls N accumulation in lake waters thereby influences the CH4 production process. The microbial community analyses revealed that the distribution of picocyanobacteria (i.e., Synechococcus), which produce CH4 in oxic conditions, was closely related to the vertical distribution of dissolved CH4 and SMM formation. Moreover, a cross-lake comparison showed that lakes with a more abundant Synechococcus population exhibited a greater development of the SMM, suggesting that these microorganisms are the most likely cause of methane production. Thus, we conclude that the stoichiometric balance between DOC and DIN might cause the cascading responses of biogeochemical processes, from N depletion to picocyanobacterial domination, and subsequently influence SMM formation in lake ecosystems.


Introduction
In recent decades, considerable efforts have been made to elucidate the source of methane (CH 4 ) emissions, because the atmospheric concentration of this greenhouse gas has been increasing [1][2][3].
are located in temperate and cool-temperate regions (35 • 11' N-42 • 36' N latitudes) and are mainly classified as caldera lakes, dammed lakes, or tectonic lakes (Table S1). The study lakes are deep stratified lakes, with maximum depths (31-363 m) deeper than the depth of 1% of the irradiance at the surface (i.e., a surrogate for the compensation depth) during the study period. The lake trophy was classified as mesotrophic, oligotrophic, or ultraoligotrophic conditions according to chlorophyll a and TP concentrations (Table S2, [30]). Most of the lakes are monomictic or dimictic and rarely freeze over; except Lake Hibara, where the surface water generally freezes during the December-March period. The watersheds of the study lakes are mostly covered with deciduous and coniferous forests (range = 67-99%). Agricultural lands such as rice paddies and farmlands occupy 20% and 11% of the watershed areas of Lake Toya and Lake Nojiri, respectively. Grassland vegetation accounts for 11% of the watershed area of Lake Hibara, while in Lake Ashinoko, residential areas (9%) contribute to the land cover to some extent.
At the deepest point of each study lake, we established a sampling site and measured the vertical profiles of environmental variables such as water temperature ( • C), dissolved oxygen (DO) concentration (DO, mgO 2 /L) and saturation (%), pH, and photosynthetically active radiation (PAR, µmol m −2 s −1 ) by using a multi-parameter sonde (YSI ProDSS, YSI-Nanotech, Tokyo, Japan), a DO meter (HQ40d, Hach, Loveland, CO, USA), a pH meter (Orion3-Star, Thermo Scientific, Chelmsford, MA, USA), and a submersible spherical quantum sensor (LI-193, LI-COR, Inc., Lincoln, NE, USA), respectively. For each sampling site, lake water samples were collected at 10 depths, from the surface to hypolimnetic waters, by using a 6 L Van Dorn sampler. For determination of the dissolved CH 4 concentration, the collected lake water was siphoned into two 30 mL glass vials and sealed with a butyl rubber stopper and an aluminum crimp. Then, 0.2 mL of saturated HgCl 2 was added to each vial as a preservative and stored in a dark place at room temperature until analysis. For water quality analyses, the collected lake water was filtered with a 100 µm nylon mesh to remove large particles (e.g., zooplankton), and the filtrates were stored in polypropylene bottles for total nitrogen (TN) and total phosphorus (TP) analyses. The remaining filtrate water was further filtered through pre-combusted GF/F glass fiber filters (Whatman plc., Maidstone, Kent, UK) for dissolved nutrient (i.e., SRP, NH 4 , NO 2 , and NO 3 ) and dissolved organic carbon (DOC) analyses; the former samples were stored in polypropylene bottles, while the latter samples were preserved in pre-combusted amber glass bottles. All of the water samples for nutrient and DOC analyses were kept at −20 • C until analyzed. For analyzing dissolved manganese (Mn) concentrations, the 100 µm mesh-filtered lake waters were further filtered with cellulose acetate disposable filters (pore size 0.45 µm), and nitric acid (final concentration = 0.1 M HNO 3 ) was added to prevent the precipitation and adsorption of metals. Then, the acidified water samples were stored in a refrigerator.
For determination of the vertical profile of phytoplankton biomass and planktonic microbe density, we also collected lake waters from the 10 sampling depths and performed chlorophyll a (Chl a) measurements and CARD-FISH (catalyzed reporter deposition fluorescence in situ hybridization) analyses. We collected phytoplankton pigment samples by filtering 50−500 mL of lake water with GF/F glass fiber filters, and then stored the samples at −20 • C until analysis. For the CARD-FISH analyses, an aliquot of lake water samples was fixed in a 0.2 µm-filtered paraformaldehyde solution (final concentration = 1%, v/v) and stored at 4 • C for <24 h. To analyze the community structure of planktonic eukaryotic algae and prokaryotic bacteria, we collected additional lake water samples from each epilimnion (0 m depth), metalimnion (7−15 m depth), and hypolimnion (26−60 m depth). The samples for the enumeration of eukaryotic algae were fixed with a Lugol solution (final concentrations, 0.2%−0.4%) and stored in the dark at room temperature, while the samples for bacterial community structure were collected on 0.22 µm filters (Sterivex filter cartridges, Millipore, Billerica, MA, USA) and stored at −20 • C until analysis.
To identify the presence of phosphonic compounds in microbial cells, a two-dimensional NMR ( 1 H-31 P heteronuclear multiple bond correlation, 1 H-31 P HMBC) analysis was performed for the suspended particles in lake water samples according to the previous study [31]. A large amount of lake water (68−201 L) was collected from the metalimnion of each lake with the Van Dorn sampler and filtered with a 100 µm nylon mesh to remove large particles. Then, the suspended particulate matter was collected onboard by filtering the collected lake water through a 2 µm cartridge filter (MCP-HX-C10S, Advantec, Inc., Tokyo, Japan) by using a peristaltic pump (Masterflex, L/S VFP002, Cole-Parmer, Vernon Hills, IL, USA). The filter samples were stored at −20 • C until analysis.
We also surveyed all of the major tributary streams identifiable on topographical maps for each study lake (Table S3), except for two inaccessible tributaries (Yanagisawa Stream of Lake Chuzenji and Okotanpe Stream of Lake Shikotsu). Near the downstream end of each tributary inflow, we collected river water samples to analyze the same environmental variables (water temperature, DO, pH, dissolved CH 4 , nutrients, DOC, and Mn) as those in the lake survey. In addition, the discharge (m 3 /s) of tributary streams was measured by the midsection method using a portable current meter (CR-7WP; Cosumo Riken, Inc., Kashiwara, Japan).

Chemical Analyses
Dissolved CH 4 concentrations of the collected lake waters were quantified by the headspace equilibration method using a gas chromatograph with an FID detector (GC-FID, GC-8A, Shimadzu Corp., Kyoto, Japan). The headspace phase created by 3 mL of high-purity He gas (>99.9999%) in the vial was equilibrated with the aqueous phase. The pCH 4 in the equilibrated headspace was then analyzed by GC-FID. The dissolved CH 4 concentration (nmol/L) in the water samples were determined from the amounts of gasses in both the equilibrated headspace and the aqueous phase; the latter was calculated by using pCH 4 in the headspace and Henry's law constant.
The DOC concentration of lake water samples was measured using a total organic carbon analyzer (TOC-L CPN, Shimadzu Corp., Kyoto, Japan). The ammonium concentration (NH 4 , µM) was quantified spectrophotometrically using the indophenol method [32]. Nitrate (NO 3 ) and nitrite (NO 2 ) concentrations were determined through the second-derivative UV spectrophotometric method [33] and Bendschneider and Robinson's method [34], respectively. Dissolved inorganic nitrogen (DIN) concentrations were determined as the sum of NH 4 , NO 2 , and NO 3 . The total nitrogen (TN) concentration of the lake water samples was determined by the alkaline persulfate digestion method, followed by the absorbance measurement (at 220 nm) for the resultant nitrate by UV spectrophotometry. Soluble reactive phosphorus (SRP) concentrations were measured spectrophotometrically using the molybdenum blue method. Likewise, the dissolved total phosphorus (DTP) and total phosphorus (TP) concentrations were quantified using the molybdenum blue method after persulfate digestion for the GF/F-filtered lake waters and unfiltered lake waters, respectively. For the phosphomolybdic acid measurements for SRP, DTP, and TP analyses, we increased the sulfuric acid concentration by a factor of two to minimize the effect of silicate interference [35]. The dissolved organic phosphorus (DOP) concentration was determined as the difference between DTP and SRP concentrations. Dissolved Mn concentrations were analyzed using inductively coupled plasma optical emission spectroscopy (SPS 3520 UV-DD, Hitachi High-Tech Science Corp., Tokyo, Japan).

Analyses for Phytoplankton Biomass and Community Structure
The Chl a sample pigments collected on GF/F glass fiber filters were extracted with N, N-dimethylformamide (6-mL) in the dark for 24 h. The Chl a content in lake water (µg-chl a/L) was then measured by the Welschmeyer method using a fluorometer (Trilogy, Turner Designs, San Jose, CA, USA). To identify the community structure of eukaryotic algae, we identified phytoplankton species to the lowest taxonomic level and enumerated the cell numbers of each algal species for the epilimnion, metalimnion and hypolimnion of each study lake by microscopic enumeration. We also enumerated red autofluorescent algal cells on the CARD-FISH filter (see below) under blue excitation (470 nm, Filter 09, Carl Zeiss Microscopy GmbH, Jena, Germany) to estimate algal cell density (cells/mL).

Analyses for Planktonic Bacterial Density and Community Structure
To analyze the distribution of specific microbes related to SMM formation, the CARD-FISH analysis was performed for the water samples collected from each sampling depth. The aliquot of lake water samples fixed with paraformaldehyde solution was filtered through a white polycarbonate membrane filter (type GTTP, pore size 0.2 µm, Millipore, Billerica, MA, USA) and a brown polycarbonate membrane filter (type GTBP, pore size 0.2 µm, Millipore, Billerica, MA, USA) for the enumeration of CARD-FISH-and DAPI-stained bacterial cells, respectively. We applied the oligonucleotide probes (labeled with horseradish peroxidese) specific for eubacteria (EUB338), type I MOB (Mg84 + Mg705), Synechococcus (405_Syn) and archaea (ARCH915) for the filters at 46 • C for > 8h for hybridization (Table S4) because these microbial groups have been reported to influence the dissolved methane profile in lakes [15,[36][37][38]. FITC-labeled tyramide solution and DAPI solution were then applied to the filter sections for the enumeration of cells on a fluorescence microscope (Zeiss Axio Lab.A1, Carl Zeiss Microscopy GmbH, Jena, Germany) under blue excitation (470 nm, Filter 10) and UV excitation (365 nm, Filter 01), respectively. We also enumerated the autofluorescent cyanobacterial cells under green excitation (545 nm, Filter 43) to estimate the density of cyanobacteria. At least 1000 cells were counted to estimate the cell density of the target bacteria (cells/mL).
To analyze the community structure of planktonic bacteria, 16S rRNA gene amplicon sequencing was performed for the filter samples collected at the epilimnion, metalimnion and hypolimnion of each lake. From the Sterivex filter samples, DNA was extracted using the PowerWater Sterivex DNA Isolation Kit (Qiagen, Hilden, Germany). From the DNA samples, sequencing libraries were prepared according to the Illumina 16S Metagenomic Sequencing Library Preparation protocol and then subjected to Illumina MiSeq sequencing. The resulting reads were processed with QIIME, an open-source bioinformatics platform for microbial community analysis [39], as follows. Quality-filtered sequences were clustered using the UCLUST algorithm with an identity threshold of 97%, to generate operational taxonomic units (OTUs). The representative sequences of the OTUs were subjected to chimera-checking to eliminate chimeric OTUs. The representatives were then used for taxonomic assignment of the OTUs via the Greengenes database, version 13_8. Greengenes annonation includes eukaryotic organelles (i.e., chloroplasts and mitochondria) as a taxonomic group in bacterial phyla. As we focused on bacterioplankton communities in the 16S rRNA gene analyses, we removed the chloroplast and mitochondria reads from the amplicon library for data analyses. The nucleotide sequences of the 16S rRNA gene obtained in the present study are available in the NCBI/SRA with the accession number of PRJNA599317.

1 H-31 P NMR Analyses for Suspended Particles
For NMR analysis, we extracted the organic phosphorus fraction from the suspended particles collected on the 2 µm cartridge filter [31]. We separated the pleated filter paper from the polypropylene filter cartridge by using a cutter and tweezers. The separated filter papers were then placed into a Teflon vial with 50 mL Milli-Q water, and the vials were subsequently immersed for approximately 1 h in a water bath set at 60 • C for the extraction [40]. The extracts were then filtered through a GF/F glass fiber filter followed by membrane filter filtration (PVDF: nominal pore size 0.45 µm, Millex, Millipore, Billerica, MA, USA). The filtrates were frozen (−20 • C) and lyophilized to concentrate the extracts for NMR analysis. Just before the NMR analyses, the lyophilized extracts were re-dissolved in 0.9 mL of EDTA solution (20 mmol L −1 ) with D 2 O (0.1 mL, Fujifilm Wako Pure Chemical Corp., Osaka, Japan) and shaken at room temperature for 10 minutes. The solution was then filtered through a membrane filter (nominal pore size: 0.45 µm, Millex, Millipore, Billerica, MA, USA) and transferred into a 5 mm NMR tube. The 1 H-31 P NMR spectrum of the extracts was recorded at 500.2MHz (for 1 H) and 202.5 MHz (for 31 P) via a JNM-ECA 500 spectrometer (JEOL Ltd., Tokyo, Japan), equipped with a 5 mm auto-tune probe. These spectra were externally referenced to 4,4-dimethyl-4-silapentane-1-sulfonic acid (DSS) and D 3 PO 4 (δ = 0 ppm), respectively. The datasets of 1 H-31 P HMBC spectra were acquired with 2048 and 256 data points and sweep widths of 9384 and 20,259 Hz for the x and y dimensions, respectively. The obtained spectra were processed with a 2 Hz line broadening using Delta, version 5.

Data Analyses
A generalized linear model (GLM) was developed to identify environmental variables that affect the dissolved CH 4 concentrations in the study lakes. We included physical (depth, water temperature, and percentage of surface irradiance), chemical (DIN, SRP, DOC, DOP, and TN/TP), and biological (Chl a) variables as the explanatory variables. Two variables (i.e., SRP and Chl a) were transformed as log 10 (x + 1) because the data were not normally distributed by visual examination. A GLM analysis with a gamma error distribution and an inverse link function was used to analyze the factors that influenced the response variable (dissolved CH 4 ). We built candidate models containing all possible combinations of the nine explanatory variables (2 9 − 1 = 511 models) and selected the best model based on the AIC value. The relative importance of an explanatory variable was evaluated by the sum of Akaike weights (w) over all models, including that variable [41]. To identify the indirect causal relationships of the environmental variables with SMM formation, we also constructed GLMs for the environmental variables that most strongly affected the dissolved CH 4 concentrations. When we constructed the GLM for DOP values, we transformed the response variable as DOP + 1 because the gamma regression analysis did not allow for the inclusion of variables with negative or zero values. The GLM analyses and model-selection procedures were performed using the MuMIn package, version 1.40.0 [42] and the MASS package, version 7.3-48 [43] of the statistical R software 3.3.3 (R Development Core Team 2017).
For the community analyses of planktonic bacteria, 16S rRNA gene sequence read counts were used to calculate the relative abundance of each OTU for eplimnetic, metalimnetic, and hypolimnetic bacterial communities, by assuming that the amplicon reads could be used as a surrogate for the abundance of each bacterial population (total number of amplicon reads = 8848-29,331 reads per sample). In this calculation, we removed the taxonomically unassigned OTUs from the analyses, as the percentage of their amplicon reads accounted only for <1.6% of the total reads. We also confirmed a significant positive correlation between the total number of bacterial read counts per unit of lake volume and the cell density of eubacteria obtained from the CARD-FISH analyses (EUB338, r = 0.50, p = 0.015). Finally, the taxonomic composition data were combined into the phylum level for the analyses. Non-metric multidimensional scaling (NMDS) analysis was performed to identify the relationships between bacterioplankton community composition and the environmental variables of the epilimnion, metalimnion, and hypolimnion of each study lake. To obtain an NMDS ordination with a low stress value at <0.2 (an indication of well-represented data by a two-dimensional representation [44]), the phyla with a relative abundance of <0.4% of the total read counts when all of the samples were combined were grouped into a category as "Others". For the environmental variables, the average values of 10 variables (i.e., depth, temperature, percentage of surface irradiance, CH 4 , Chl a, DIN, SRP, DOC, DOP, TN/TP) were determined for each of epilimnion, metalimnion and hypolimnion, and were then incorporated in the analyses. The NMDS analysis (based on the Bray-Curtis dissimilarity) was performed by using the vegan 2.4-5 package [45] of the R software 3.3.3.
Finally, we examined the relationship between the degree of SMM development and lake environmental variables to identify the factors controlling the oxic CH 4 peak. According to our previous study [15], the SMM is herein defined as the layer of the local CH 4 maximum that forms below the surface water. The degree of SMM development for each study lake was quantified by peak SMM, namely the maximum CH 4 concentration in the SMM minus the atmospheric equilibrium CH 4 concentration at that depth. The peak SMM for the study lakes was regressed against the environmental variables responsible for SMM formation.

SMM Formation in Lakes
The vertical profiles of dissolved CH 4 concentrations revealed that the SMM developed in aerobic layers in all of the nine study lakes ( Figure 1). Further, peak SMM concentrations (range = 62.5-592.6 nM) showed the highest values in the vertical CH 4 profiles for all of the study lakes. The SMM peak tended to occur within the metalimnion, or in depths adjacent to the metalimnion where DO saturation levels were high (range = 85-131%). In addition, the peak depth and vertical profile of dissolved CH 4 (i.e., subsurface peak of CH 4 ) were similar to those of DO saturation, except in Lake Nojiri.
Water 2019, 11, x FOR PEER REVIEW 7 of 19 SMM peak tended to occur within the metalimnion, or in depths adjacent to the metalimnion where DO saturation levels were high (range = 85%-131%). In addition, the peak depth and vertical profile of dissolved CH4 (i.e., subsurface peak of CH4) were similar to those of DO saturation, except in Lake Nojiri. Figure 1. Vertical profiles of dissolved CH4 concentrations in the nine study lakes of Japan during the July to September in 2016-2017 period. AS: Lake Ashinoko, MO: Lake Motosu, SA: Lake Saiko, AO: Lake Aoki, NO: Lake Nojiri, CH: Lake Chuzenji, HI: Lake Hibara, TO: Lake Toya, SH: Lake Shikotsu. Red and gray circles refer to the dissolved methane (nM) and oxygen (%) concentrations, respectively. Dotted lines denote the atmospheric equilibrium concentration of dissolved CH4 in lake water (3.7-7.2 nM). Shaded areas denote the range of the thermocline.
Regarding the other environmental variables, water temperature and percentage of surface irradiance were highest at the lake surface and decreased with depth; this pattern did not correspond with the vertical profile of the dissolved CH4 concentration for any of the study lakes ( Figures S2 and  S3). Although the overall pattern of Chl a and DOP profiles did not match the CH4 profile, their peaks were located at depths similar to the CH4 peaks ( Figures S4 and S5). In most of the study lakes, DIN concentrations tended to be higher in the hypolimniion and lower in both the epilimnion and metalimnion, the latter of which exhibited high CH4 concentrations ( Figure S6). Although SRP did not show a clear overall trend, the peak depth and vertical profile of SRP were similar to those of CH4 in Lake Motosu, Lake Saiko, Lake Toya and Lake Shikotsu ( Figure S7). Similarly, the vertical patterns of DOC were variable, but its concentrations tended to be higher in the metalimnion or nearby layers (except in Lake Nojiri), as observed in the CH4 profile ( Figure S8). The observed values of lake TN/TP ratio (median = 151, range = 32-875) were high relative to the Redfield ratio (N:P = 16:1) for all of the study lakes ( Figure S9). The vertical TN/TP profile tended to show a pattern opposite to that of the CH4 profile-except in Lake Shikotsu, where a close association was observed for both variables.
The cross-lake comparison showed that the development of SMM varied among the study lakes ( Figure 1). The highest value of peak SMM was observed in Lake Saiko (an oligotrophic lake), while small CH4 peaks were observed in three ultraoligotrophic lakes (Lake Motosu, Lake Toya, and Lake Shikotsu). AS: Lake Ashinoko, MO: Lake Motosu, SA: Lake Saiko, AO: Lake Aoki, NO: Lake Nojiri, CH: Lake Chuzenji, HI: Lake Hibara, TO: Lake Toya, SH: Lake Shikotsu. Red and gray circles refer to the dissolved methane (nM) and oxygen (%) concentrations, respectively. Dotted lines denote the atmospheric equilibrium concentration of dissolved CH 4 in lake water (3.7-7.2 nM). Shaded areas denote the range of the thermocline.

Relationships between Dissolved CH4 and Physicochemical Variables
Regarding the other environmental variables, water temperature and percentage of surface irradiance were highest at the lake surface and decreased with depth; this pattern did not correspond with the vertical profile of the dissolved CH 4 concentration for any of the study lakes ( Figures S2 and  S3). Although the overall pattern of Chl a and DOP profiles did not match the CH 4 profile, their peaks were located at depths similar to the CH 4 peaks ( Figures S4 and S5). In most of the study lakes, DIN concentrations tended to be higher in the hypolimniion and lower in both the epilimnion and metalimnion, the latter of which exhibited high CH 4 concentrations ( Figure S6). Although SRP did not show a clear overall trend, the peak depth and vertical profile of SRP were similar to those of CH 4 in Lake Motosu, Lake Saiko, Lake Toya and Lake Shikotsu ( Figure S7). Similarly, the vertical patterns of DOC were variable, but its concentrations tended to be higher in the metalimnion or nearby layers (except in Lake Nojiri), as observed in the CH 4 profile ( Figure S8). The observed values of lake TN/TP ratio (median = 151, range = 32-875) were high relative to the Redfield ratio (N:P = 16:1) for all of the study lakes ( Figure S9). The vertical TN/TP profile tended to show a pattern opposite to that of the CH 4 profile-except in Lake Shikotsu, where a close association was observed for both variables. The cross-lake comparison showed that the development of SMM varied among the study lakes ( Figure 1). The highest value of peak SMM was observed in Lake Saiko (an oligotrophic lake), while small CH 4 peaks were observed in three ultraoligotrophic lakes (Lake Motosu, Lake Toya, and Lake Shikotsu).

Relationships between Dissolved CH 4 and Physicochemical Variables
GLM analyses were performed to identify the physicochemical variables affecting the water-column CH 4 concentrations in lakes. In the analyses, six environmental variables (i.e., depth, water temperature, percent of surface irradiance, DIN, DOP, and TN/TP ratio) were retained in the best model selected by AIC ( Figure 2). The sum of Akaike weights (w) further identified that depth, DIN, DOP, and TN/TP ratio, were the most important variables explaining the variability of dissolved CH 4 concentrations in the lakes (Table S5). GLM analyses were performed to identify the physicochemical variables affecting the watercolumn CH4 concentrations in lakes. In the analyses, six environmental variables (i.e., depth, water temperature, percent of surface irradiance, DIN, DOP, and TN/TP ratio) were retained in the best model selected by AIC (Figure 2). The sum of Akaike weights (w) further identified that depth, DIN, DOP, and TN/TP ratio, were the most important variables explaining the variability of dissolved CH4 concentrations in the lakes (Table S5).
We observed a strong negative nonlinear relationship between DIN and dissolved CH4 concentrations when we combined the data from all nine study lakes ( Figure 2 and Table S6). The CH4 concentration increased rapidly with a decrease in the DIN concentration. Similar curvilinear patterns were observed for depth and the TN/TP ratio; the lake waters with shallower depths and a lower TN/TP ratio increased the water-column CH4 concentrations. On the other hand, a positive effect of DOP on dissolved CH4 concentrations was observed for the study lakes.
To explore the potential factors controlling the chemical variables that strongly affect dissolved CH4, we further constructed GLMs for DIN, DOP, and TN/TP (Table S5). The results showed that DOC concentration was the best explanatory variable for all three influential variables ( Figure 3 and Table S6). DOC exhibited negative nonlinear relationships with both DIN and the TN/TP ratio. In contrast, DOC showed a positive relationship with DOP concentrations. Relationships between environmental variables and dissolved CH 4 concentrations in the study lakes. Univariate regression lines are shown for the explanatory variables that were retained in the best generalized linear models (GLMs) with a gamma distribution of the response variable. : Lake Ashinoko, : Lake Motosu, : Lake Saiko, : Lake Aoki, : Lake Nojiri, : Lake Chuzenji, : Lake Hibara, : Lake Toya, : Lake Shikotsu.
We observed a strong negative nonlinear relationship between DIN and dissolved CH 4 concentrations when we combined the data from all nine study lakes ( Figure 2 and Table S6). The CH 4 concentration increased rapidly with a decrease in the DIN concentration. Similar curvilinear patterns were observed for depth and the TN/TP ratio; the lake waters with shallower depths and a lower TN/TP ratio increased the water-column CH 4 concentrations. On the other hand, a positive effect of DOP on dissolved CH 4 concentrations was observed for the study lakes.
To explore the potential factors controlling the chemical variables that strongly affect dissolved CH 4 , we further constructed GLMs for DIN, DOP, and TN/TP (Table S5). The results showed that DOC concentration was the best explanatory variable for all three influential variables ( Figure 3 and Table S6). DOC exhibited negative nonlinear relationships with both DIN and the TN/TP ratio. In contrast, DOC showed a positive relationship with DOP concentrations.

Community Analyses of Planktonic Algae and Bacteria
Microscopic enumeration for eukaryotic algae identified six phytoplankton classes present in the study lakes (Table S7). Among these phytoplanktonic algal groups, Chlorophyceae or Bacillariophyceae predominated in the epilimnion of most of the study lakes except for Lake Toya, where Dinophyceae (53.6%) was the dominant group. Similarly, Chlorophyceae and/or Bacillariophyceae tended to predominate in the algal communities in the metalimnion, although Cryptophyceae and Dinophyceae were relatively abundant in the metalimnion of Lake Ashinoko. In the hypolimnion, Bacillariophyceae accounted for 54%-97% of the relative abundance of algal communities, although Chlorophyceae was relatively abundant in the hypolimnion community of Lake Chuzenji (45.8%). There was no phytoplanktonic algal group that consistently exhibited a metalimnetic peak in terms of relative abundance for any of the study lakes.

Community Analyses of Planktonic Algae and Bacteria
Microscopic enumeration for eukaryotic algae identified six phytoplankton classes present in the study lakes (Table S7). Among these phytoplanktonic algal groups, Chlorophyceae or Bacillariophyceae predominated in the epilimnion of most of the study lakes except for Lake Toya, where Dinophyceae (53.6%) was the dominant group. Similarly, Chlorophyceae and/or Bacillariophyceae tended to predominate in the algal communities in the metalimnion, although Cryptophyceae and Dinophyceae were relatively abundant in the metalimnion of Lake Ashinoko. In the hypolimnion, Bacillariophyceae accounted for 54-97% of the relative abundance of algal communities, although Chlorophyceae was relatively abundant in the hypolimnion community of Lake Chuzenji (45.8%). There was no phytoplanktonic algal group that consistently exhibited a metalimnetic peak in terms of relative abundance for any of the study lakes.
The NMDS analyses showed that the community composition of epilimnetic and metalimnetic microbes differed from that of the hypolimnetic community ( Figure 4). The results revealed that the relative abundance of Proteobacteria, Verrucomicrobia, Actinobacteria, Bacteroidetes, and Cyanobacteria tended to be high in the epilimnion and metalimnion, whereas Planctomycetes, Acidobacteria, Chloroflexi, and Chlorobi were relatively abundant in the hypolimnion. We further observed the association of community structure and environmental variables. In epilimnetic and metalimnetic layers with high CH 4 and DOP concentrations, Cyanobacteria, and Verrucomicrobia tended to be relatively abundant.
The results of epifluorescent microscopic analyses (CARD-FISH, DAPI, and autofluorescent enumeration) showed that both Cyanobacteria and Synechococcus cell density showed significant positive relationships with water-column dissolved CH 4 concentrations ( Figure 5). Synechococcus accounted for 21-96%, 64-96%, and 22-75% of the total cyanobacterial density in the epilimnion, metalimnion, and hypolimnion of the study lakes, respectively. In contrast to the patterns of Cyanobacteria and Synechococcus, phytoplankton biomass (Chl a), as well as the cell density of eubacteria (EUB338), type I MOB (Mg84 + Mg705), archaea (ARCH915), and algae, showed no significant relationships with dissolved CH 4 concentrations. metalimnetic peak in terms of relative abundance for any of the study lakes.
The results of epifluorescent microscopic analyses (CARD-FISH, DAPI, and autofluorescent enumeration) showed that both Cyanobacteria and Synechococcus cell density showed significant positive relationships with water-column dissolved CH4 concentrations ( Figure 5). Synechococcus accounted for 21%-96%, 64%-96%, and 22%-75% of the total cyanobacterial density in the epilimnion, metalimnion, and hypolimnion of the study lakes, respectively. In contrast to the patterns of Cyanobacteria and Synechococcus, phytoplankton biomass (Chl a), as well as the cell density of eubacteria (EUB338), type I MOB (Mg84 + Mg705), archaea (ARCH915), and algae, showed no significant relationships with dissolved CH4 concentrations.  Table S4). GLM with a Gaussian error distribution and an identity link function was used for the variables, while regression lines shown when significant (p < 0.05). See the legend for Figure 2 for a description of the symbols.
We then analyzed the vertical profile of Synechococcus density in relation to the CH4 profiles, which showed that both profiles were closely associated with each other (Figure 6). In addition, a negative curvilinear pattern was found for the relationship of Synechococcus cell density with both depth and DIN concentration, as was observed in the patterns of dissolved CH4 concentrations ( Figure S10). In contrast, a positive nonlinear relationship between Synechococcus cell density and DOP concentrations was found.  Table S4). GLM with a Gaussian error distribution and an identity link function was used for the variables, while regression lines shown when significant (p < 0.05). See the legend for Figure 2 for a description of the symbols.
We then analyzed the vertical profile of Synechococcus density in relation to the CH 4 profiles, which showed that both profiles were closely associated with each other (Figure 6). In addition, a negative curvilinear pattern was found for the relationship of Synechococcus cell density with both depth and DIN concentration, as was observed in the patterns of dissolved CH 4 concentrations ( Figure S10). In contrast, a positive nonlinear relationship between Synechococcus cell density and DOP concentrations was found.

Relationship between Peak SMM and Environmental Variables
A correlation analysis between the peak SMM and environmental variables revealed that the peak SMM was positively and linearly correlated with lake DOC, DOP, and Synechococcus cell density observed at the same depth ( Figure 7 and Table S6). Lakes with higher concentrations of both DOC and DOP and higher Synechococcus cell density tended to develop a larger SMM peak. Moreover, a negative nonlinear relationship of peak SMM was observed for both DIN concentration and the TN/TP ratio, implying that lakes with lower DIN and TN/TP values experienced a rapid increase in peak SMM.

Relationship between Peak SMM and Environmental Variables
A correlation analysis between the peak SMM and environmental variables revealed that the peak SMM was positively and linearly correlated with lake DOC, DOP, and Synechococcus cell density observed at the same depth ( Figure 7 and Table S6). Lakes with higher concentrations of both DOC and DOP and higher Synechococcus cell density tended to develop a larger SMM peak. Moreover, a negative nonlinear relationship of peak SMM was observed for both DIN concentration and the TN/TP ratio, implying that lakes with lower DIN and TN/TP values experienced a rapid increase in peak SMM.

Relationship between Peak SMM and Environmental Variables
A correlation analysis between the peak SMM and environmental variables revealed that the peak SMM was positively and linearly correlated with lake DOC, DOP, and Synechococcus cell density observed at the same depth (Figure 7 and Table S6). Lakes with higher concentrations of both DOC and DOP and higher Synechococcus cell density tended to develop a larger SMM peak. Moreover, a negative nonlinear relationship of peak SMM was observed for both DIN concentration and the TN/TP ratio, implying that lakes with lower DIN and TN/TP values experienced a rapid increase in peak SMM.  Figure 7. Relationships between lake environmental variables and the peak subsurface methane maximum (SMM) for the study lakes. For the relationships of both dissolved inorganic nitrogen (DIN) and TN/TP, the regression lines were estimated using generalized linear models (GLMs) with a gamma error distribution and an inverse link function for each variable, while the GLM with a Gaussian error distribution and an identity link function was used for the other variables. For explanatory variables, the values at the peak SMM depths were used for the analyses.

Identification of Phosphonates in Suspended Particles
Particulate suspended matter collected from the metalimnion was characterized by a two-dimensional 1 H-31 P NMR analysis. The spectral region of 20-31 ppm in 31 P NMR analysis is known to represent the chemical shifts of phosphonate compounds [25,46]. Our mass spectra at 1.2 ppm ( 1 H) and 25 ppm ( 31 P) revealed the presence of phosphonate in the suspended particles of the metalimnion in Lake Toya (Figure 8).
Water 2019, 11, x FOR PEER REVIEW 12 of 19 and TN/TP, the regression lines were estimated using generalized linear models (GLMs) with a gamma error distribution and an inverse link function for each variable, while the GLM with a Gaussian error distribution and an identity link function was used for the other variables. For explanatory variables, the values at the peak SMM depths were used for the analyses.

Identification of Phosphonates in Suspended Particles
Particulate suspended matter collected from the metalimnion was characterized by a twodimensional 1 H-31 P NMR analysis. The spectral region of 20-31 ppm in 31 P NMR analysis is known to represent the chemical shifts of phosphonate compounds [25,46]. Our mass spectra at 1.2 ppm ( 1 H) and 25 ppm ( 31 P) revealed the presence of phosphonate in the suspended particles of the metalimnion in Lake Toya (Figure 8).

Discussion
In the present study, we observed the formation of the CH4 maximum in oxic subsurface layers (i.e., SMM) within or near the metalimnion in all of the study lakes (Figure 1). The seasonal development of metalimnetic CH4 peaks during the stratification period was also reported by previous studies in deep freshwater lakes [8][9][10]15] as well as in pelagic marine ecosystems [16,17]. These results suggest that the development of the CH4 peak in aerobic subsurface waters may be a common phenomenon in aquatic ecosystems.

Discussion
In the present study, we observed the formation of the CH 4 maximum in oxic subsurface layers (i.e., SMM) within or near the metalimnion in all of the study lakes (Figure 1). The seasonal development of metalimnetic CH 4 peaks during the stratification period was also reported by previous studies in deep freshwater lakes [8][9][10]15] as well as in pelagic marine ecosystems [16,17]. These results suggest that the development of the CH 4 peak in aerobic subsurface waters may be a common phenomenon in aquatic ecosystems.
SMM formation was unlikely to couple with the dissolution of atmospheric CH 4 because the observed CH 4 peak concentrations (67-597 nM, Figure 1) were one or two orders of magnitude higher than the dissolved CH 4 equilibrated with the atmosphere (range = 3.4-7.4 nM). Tributary inflow is known to contribute to SMM development in lakes [9]. However, we observed apparent SMM development even in lakes with no tributary streams (i.e., Lake Motosu and Lake Aoki). For the other lakes, moreover, the discharge-weighted average of the dissolved CH 4 concentrations of tributary streams (range = 8.8-88 nM) was generally lower than the peak CH 4 concentrations-except in Lake Ashinoko, where there is a small tributary (discharge < 0.02 m 3 /s) with a high CH 4 concentration (Table S3). Another potential source of SMM formation might be the transportation of CH 4 from anoxic littoral and profundal sediments [20,21]. However, the vertical profile of dissolved Mn concentrations, a tracer of water mass transported from anoxic environments [9,15], showed a pattern unrelated to the CH 4 profiles ( Figure S11). Moreover, the vertical gradient of CH 4 profiles suggests no substantial contribution of hypolimnetic CH 4 to the SMM of any of the study lakes ( Figure 1). Therefore, we argue that the diffusion of anoxic CH 4 from anoxic littoral and profundal sediments was also unlikely to be the source of the CH 4 supersaturation observed in our study.
Recent studies have suggested that the microbial degradation of phosphonates (such as MPn and 2-AEP) can explain aerobic methane production in oligotrophic waters [12,13,15,16,24,25]. The marine and freshwater microbial community (e.g., Trichodesmium, Pseudomonas, Synechococcus, SAR11) can express the C-P lyase operon (phn genes) to utilize phosphonic compounds under phosphate-starved conditions [15,16,24,25]. As a result of the cleavage of C-P bonds of phosphonates, these microbes can produce CH 4 . In fact, the batch-culture experiments of lake water amended with MPn confirmed in situ aerobic CH 4 production by planktonic bacteria in one of the study lakes (Lake Saiko) [15]. Considering such circumstantial evidences, as well as the fact that the DO saturation profile resembled the dissolved CH 4 pattern (Figure 1), we argue that photosynthesis-related biogeochemical processes mediated by planktonic microbes may be relevant to the development of SMM.
The relationships between dissolved CH 4 and limnological variables suggest that DOC controls the availability of DIN, DOP, and the TN/TP ratio, thereby influencing the water-column CH 4 concentration in the study lakes (Figures 2 and 3). In particular, the remarkable negative nonlinear relationships of DOC with DIN and the subsequent nonlinear effects of DIN on dissolved CH 4 concentrations were found. The nonlinear relationship between DOC and nitrate (NO 3 ) has recently been reported from a wide variety of aquatic environments along a hydrologic continuum from soils to streams and lakes to coastal and pelagic ocean ecosystems [47][48][49][50]. It has been hypothesized that NO 3 accumulation may occur in aquatic environments with low DOC concentrations due to an organic C limitation for assimilation and/or denitrification by heterotrophic microbes [50]. In contrast, NO 3 depletion may occur in aquatic environments with high DOC concentrations due to the sufficient supply of organic C for heterotrophs. These patterns imply that stoichiometric controls of DOC-NO 3 relationships over microbial activity may regulate the fate of N in aquatic ecosystems. In fact, NO 3 is the dominant species of DIN in the study lakes (median = 85%, range = 0-98%). Moreover, even when we separately analyzed individual lakes, a significant negative correlation of DOC and DIN was found for each of them (r = −0.69-−0.90, p < 0.05), except Lake Shikotsu (r = −0.44, p > 0.05). Furthermore, a strong nonlinear effect of DOC on the TN/TP ratio was also evident (Figure 2), suggesting that the availability of organic carbon controls N accumulation in the water column, thereby influencing the elemental stoichiometric balance of nitrogen and phosphorus in lakes.
Previous studies explained that aerobic CH 4 production from phosphonate decomposition by planktonic microbes occurs under P-starved environments because the organisms utilize phosphonates as an alternative P source [12,13,15,16,24]. However, the present study showed that higher CH 4 concentrations tended to occur in lower DIN concentrations; no clear pattern was found for the relationship between SRP and CH 4 ( Figure 2). These patterns seemed to be intuitively paradoxical, as a lower N availability was likely to stimulate aerobic CH 4 production. However, the observed TN/TP ratio in the study lakes (median = 151, range = 32-875, Table S2) was much higher than the Redfield ratio (N/P = 16), suggesting that lake productivity may be more P-limited than N-limited in these lakes. Therefore, we predict that P-starved planktonic microbes utilize phosphonate compounds to form SMM in the oxic waters of the study lakes. The positive relationship between DOP, which may contain dissolved organic phosphonate compounds, and CH 4 concentrations, also supports this argument ( Figure 2). Moreover, the NMR analysis revealed the presence of particulate phosphonates in the metalimnion of Lake Toya, where SMM formation was observed (Figure 8). A previous study also detected 2-aminoethylphosphonate (2-AEP), a possible precursor of aerobic methane production [51], from the suspended particles in the epilimnion of Lake Saiko [31]. Therefore, we hypothesized that planktonic bacteria in deep unproductive lakes may store P as phosphonate molecules and utilize such intracellular P compounds via the enzymatic liberation of C-P bonds under P-limited environments [15]. Although we were unable to detect phosphonates in the suspended fractions from the other study lakes, this is probably due to the insensitivity of NMR; prolonged data acquisition and/or a larger volume of lake water filtration may be necessary to obtain the spectra. In addition, further development of quantitative analyses for the dissolved species of phosphonate compounds [52] is necessary to understand the phosphonate and CH 4 dynamics in lakes.
The 16S rRNA gene amplicon sequencing analyses and the CARD-FISH analyses revealed that cyanobacteria, especially Synechococcus, were related to dissolved CH 4 concentrations and SMM formation. The composition of cyanobacteria estimated from amplicon reads was associated with the gradient of the CH 4 concentration (Figure 4). Moreover, a close association of the vertical CH 4 profile with Synechococcus distribution was observed across the lakes (Figure 6), as found in a previous study [15]. The batch-culture experiments conducted in previous studies also showed that Synechococcus strains have the ability to produce CH 4 in oxic conditions as a result of phosphonate (such as MPn and 2-AEP) utilization [15,51] or unknown photosynthesis-related processes [26]. Although other heterotrophic bacteria, such as Limnohabotans and Pseudomonas, are also known to carry the C-P lyase (phnJ) gene [13,15], the strong correlation between the CH 4 concentration and Synechococcus cell density ( Figure 5), as well as the association of the DO saturation profile with the CH 4 peak (Figure 1), suggested that Synechococcus (photosynthetic cyanobacteria) may be one of the potential drivers of SMM formation in deep freshwater lakes.
Synechococcus, a group of ubiquitous freshwater picocyanobacteria, often predominate in the planktonic communities in the epilimnion and metalimnion of unproductive lakes during the summer stratification period, where nutrient depletion occurs due to the prevention of vertical water mixing and nutrient exhaustion by competitors, such as eukaryotic algae [53][54][55]. Our previous studies also indicated that the development of nutrient-depleted environments by stratification during midsummer (i.e., low N and P conditions) may favor the growth of Synechococcus because their small cell size is advantageous for nutrient uptake under oligotrophic conditions due to efficient nutrient diffusion per unit of cell volume [15,55]. The present results add a further potential scenario for picocyanobacterial dominance in the metalimnion; that is, the reduced nutrient availability associated with the increase in DOC loading in lakes may promote the population growth of Synechococcus. This may be the reason for the cascading negative responses of biogeochemical elements from DOC to DIN, and from DIN to both Synechococcus and CH 4 , that were observed in our study. In fact, a cross-lake comparison revealed that lakes with relatively high DOC and low DIN concentrations enhanced the degree of SMM development (Figure 7). Likewise, lakes with more abundant Synechococcus populations and a higher DOP pool promoted the development of SMM (Figure 7). Therefore, we conclude that the stoichiometric balance between DOC and DIN in the water column may regulate nutrient availability for bacterioplankton communities (e.g., Synechococcus) and may subsequently control CH 4 production and SMM formation in deep freshwater lakes (Figure 9). Although the variability of DOC and DIN concentrations could also be explained by the extracellular release of DOC and DIN uptake by phytoplankters [28,29], the relationship between phytoplankton biomass (Chl a) and DOC concentration, as well as that between Chl a and DIN, was unclear ( Figure S12). Moreover, the vertical profiles of DOC and Synechococcus density were not closely associated with each other for the study lakes ( Figure S13), except Lake Saiko and Lake Hibara. Therefore, we argue that phytoplankton (including Synechococcus) per se may not be a major factor controlling the DOC and DIN variability, even though further studies are necessary to clarify the causal relationships between microbial activities and lake biogeochemical conditions. supply of nutrients from hypolimnion [56], which might result in a low DIN/DOC ratio in shallow water, as observed in the present study ( Figure 9). If this is the case, the projected alteration in the lake mixing regime from dimictic to monomictic, or monomictic to meromictic, under climate change might increase the aerobic CH4 production and emission from the nutrient-depleted oxic lake waters. Further studies are required to precisely predict the effect of lake-mixing regime on SMM formation. Figure 9. Effects of the stoichiometric balance between DIN and DOC on nitrogen availability, Synechococcus density and dissolved CH4 concentration in the study lakes. See the legend for Figure 2 for a description of the symbols.

Conclusion
The formation of SMM might have a significant impact on methane emissions from deep freshwater lakes because methanotrophs generally consume the majority of CH4 and prevent fugitive methane emissions from deep profundal sediments or hypolimnetic zones in lakes [36,37]. In contrast, CH4, produced by aerobic methanogenesis within the SMM of the metalimnetic layer, might easily leak into the atmosphere because of the proximity of CH4 production sites to the lake surface. Therefore, the present study predicted that DOC loading triggers the cascading responses of biogeochemical processes, from N depletion to picocyanobacterial domination, which may promote CH4 production and emission to the atmosphere from the SMM layer. Although previous studies also reported that DOC might influence the CH4 emissions of freshwater ecosystems [57,58], these studies considered DOC as a direct substrate for anaerobic methanogenesis in sediments. In contrast, our study suggested that organic C exerts stoichiometric control over CH4 production in the oxic layer of lakes. DOC concentrations are now expected to increase due to ongoing climate change, which may result in significant changes to the structure and functioning of lake ecosystems [48]. For example, increasing concentrations of DOC often change the light and thermal environments in water columns, thereby affecting the primary and secondary productivity of lake food webs, as well as the thermal stratification regime in lakes. An increased DOC supply can also provide energy subsidies to heterotrophic consumers and consequently influence the metabolic balance of lake ecosystems. Moreover, DOC can change the chemical environments of lake ecosystems, such as their pH and dissolved iron concentrations. We argue that, in addition to these phenomena, increased organic carbon loading under changing environments may promote aerobic CH4 production and emission from the oxic layer of deep freshwater lakes.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, supplementary data accompany this paper with Figures S1: Geographic location of nine study lakes in Japan, Figures S2: Vertical profiles of dissolved CH4 concentration (nM, red circles) and water temperature (°C, gray circles) in nine study lakes during the period from July to September in 2016-2017. AS: Lake Ashinoko, MO: Lake Motosu, SA: Lake Saiko, AO: Lake Aoki, NO: Lake Nojiri, CH: Lake Chuzenji, HI: Lake Hibara, TO: Lake Toya, SH: Lake Shikotsu. Although the present studies are conducted in monomictic or dimictic lakes where stratification occurs seasonally, enhanced SMM development might be predicted in meromictic lakes. This is because the lack of vertical mixing in such permanently stratified lakes may reduce the up-welling supply of nutrients from hypolimnion [56], which might result in a low DIN/DOC ratio in shallow water, as observed in the present study ( Figure 9). If this is the case, the projected alteration in the lake mixing regime from dimictic to monomictic, or monomictic to meromictic, under climate change might increase the aerobic CH 4 production and emission from the nutrient-depleted oxic lake waters. Further studies are required to precisely predict the effect of lake-mixing regime on SMM formation.

Conclusions
The formation of SMM might have a significant impact on methane emissions from deep freshwater lakes because methanotrophs generally consume the majority of CH 4 and prevent fugitive methane emissions from deep profundal sediments or hypolimnetic zones in lakes [36,37]. In contrast, CH 4 , produced by aerobic methanogenesis within the SMM of the metalimnetic layer, might easily leak into the atmosphere because of the proximity of CH 4 production sites to the lake surface. Therefore, the present study predicted that DOC loading triggers the cascading responses of biogeochemical processes, from N depletion to picocyanobacterial domination, which may promote CH 4 production and emission to the atmosphere from the SMM layer. Although previous studies also reported that DOC might influence the CH 4 emissions of freshwater ecosystems [57,58], these studies considered DOC as a direct substrate for anaerobic methanogenesis in sediments. In contrast, our study suggested that organic C exerts stoichiometric control over CH 4 production in the oxic layer of lakes. DOC concentrations are now expected to increase due to ongoing climate change, which may result in significant changes to the structure and functioning of lake ecosystems [48]. For example, increasing concentrations of DOC often change the light and thermal environments in water columns, thereby affecting the primary and secondary productivity of lake food webs, as well as the thermal stratification regime in lakes. An increased DOC supply can also provide energy subsidies to heterotrophic consumers and consequently influence the metabolic balance of lake ecosystems. Moreover, DOC can change the chemical environments of lake ecosystems, such as their pH and dissolved iron concentrations. We argue that, in addition to these phenomena, increased organic carbon loading under changing environments may promote aerobic CH 4 production and emission from the oxic layer of deep freshwater lakes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/2/402/s1, supplementary data accompany this paper with Figure S1: Geographic location of nine study lakes in Japan, Figure S2: Vertical profiles of dissolved CH 4 concentration (nM, red circles) and water temperature ( • C, gray circles) in nine study lakes during the period from July to September in 2016-2017. AS: Lake Ashinoko, MO: Lake Motosu, SA: Lake Saiko, AO: Lake Aoki, NO: Lake Nojiri, CH: Lake Chuzenji, HI: Lake Hibara, TO: Lake Toya, SH: Lake Shikotsu. Shaded areas denote the range of the thermocline, Figure S3: Vertical profiles of dissolved CH 4 concentration (nM, red circles) and the percentage of surface irradiance (%, gray circles) in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S4: Vertical profiles of dissolved CH 4 (nM, red circles) and chlorophyll a (µg/L, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S5: Vertical profiles of dissolved CH 4 (nM, red circles) and DOP (µM, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S6: Vertical profiles of dissolved CH 4 (nM, red circles) and DIN (µM, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S7: Vertical profiles of dissolved CH 4 (nM, red circles) and SRP (µM, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S8: Vertical profiles of dissolved CH 4 (nM, red circles) and DOC (µM, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S9: Vertical profiles of dissolved CH 4 concentration (nM, red circles) and TN/TP (ratio, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S10: Relationships between Synechococcus cell density (cells/mL) and lake physicochemical variables in study lakes. Univariate regression lines were shown for the explanatory variables that were retained in the best model of GLMs with gamma distribution of response variable. See the legend of Figure S10 for symbol description, Figure S11: Vertical profiles of dissolved CH 4 concentration (nM, red circles) and Mn (µM, gray circles) concentrations in nine study lakes during the period from July to September in 2016-2017. See legend of Figure S2 for the abbreviation of lake names and the descriptions of shaded areas, Figure S12: Relationships between phytoplankton biomass (log 10 (Chl a), µg/L) and two variables (DIN and DOC, µM) in study lakes, Figure S13: Vertical distribution of Synechococcus density (red circles) in relation to the DOC profile (gray circles) in the study lakes. Shaded areas denote the range of the thermocline. See legend of Figure S2 for the abbreviation of lake names, Table S1: Watershed and lake morphological variables of the nine study lakes in Japan, Table S2: Water quality variables of the nine study lakes in Japan, Table S3: Water quality variables of tributary streams for study lakes, Table S4: The CARD-FISH probes used in the study, Table S5: Relative importance for limnological explanatory variables used in the generalized linear models (GLMs) with a gamma error distribution and an inverse link function for each response variable. Relative importance was evaluated by the sum of Akaike weights (w) for each model (CH 4 , DIN, DOP and TN/TP), Table S6: Results of the generalized linear models (GLM) assessing the univariate relationship between environmental variables and dissolved CH 4 concentrations (as shown in Figure 2), the univariate relationship between DOC and chemical variables (Figure 3), and the univariate relationship between lake environmental variables and the peak subsurface maximum (SMM) in nine study lakes (Figure 7). GLMs were constructed with a gamma error distribution and an inverse link function for each variable (1/y = ax + b), except for the variables with an asterisk whose relationships were estimated with GLMs with a Gaussian error and an identity link function (y = ax + b), Table S7: Relative abundance (%) of phytoplanktonic algae in the epilimnion, metalimnion and hypolimnion of Table S6. Relative abundance (%) of phytoplanktonic algae in the epilimnion, metalimnion and hypolimnion of study lakes.