Vertical biogeochemical zonation according to patterns of energy metabolism. Total S/N/CH 4 metabolic abundances with seasonal variations ( A a) were predicted with 16S rRNA genes; specific function modules and chemotaxis ( A b–d) were predicted by metagenomic data annotated with KEGG and NCycDB. Sulfate reduction ( R 2 = 0.975, df = 5) and nitrification ( R 2 = 0.854, df = 6) were fitted with exponential model of plateau followed by one phase decay. Methanotrophy ( R 2 = 0.983, df = 4), methanogenesis ( R 2 = 0.999, df = 4), and nitrate reduction ( R 2 = 0.998, df = 4) were fitted piecewise with a linear model followed by exponential one-phase decay/association. Least squares fit was used for all these nonlinear models. Relative abundance of the NC10 class ( A c) serves as proxy for the intensity of N-damo, the nitrate/nitrite-dependent anaerobic methane oxidation. Black dots in plot A c represent the median value of the seasonal samples labeled by gray dots. OATZ oxic–anoxic transition zone, NATZ nitrate–ammonium transition zone, SMTZ sulfate–methane transition zone, NMTZ nitrate/nitrite–methane transition zone, MGZ methanogenetic zone.