In May 2014, two parallel sediment cores (CQ1: 37 cm and CQ2: 30 cm) through the sediment–water interface, were retrieved from the deepest site of the lake at a water depth of ~26 m using a gravity corer (Figure 1C; Chai et al., 2018; Zhang et al., 2022). In the field, both sediment cores were sub-sectioned by extruding vertically at every 0.5-cm interval immediately and the sedimentary lithology is mainly dominated by humus black clay. The sediment samples were transported to the laboratory and stored in the freezer before freeze-drying. The CQ2 is used for dating analysis and the CQ1 for proxy analysis including sedimentary pigment, total organic carbon (TOC), TN, carbon isotope (δ¹³C), nitrogen isotope (δ¹⁵N), and TP. The chronology is determined by ²¹⁰Pb and ¹³⁷CS activities analysis using a gamma spectrometer (Hyperpure Ge detector) at 46.5 KeV and 661.6 KeV, respectively. Chronology of CQ2 were established based on ²¹⁰Pb and ²²⁶Ra activities using the constant rate of supply (CRS) model and the chronology of core CQ1 was obtained from core CQ2 through the correlation of loss-on-ignition profiles (Chai et al., 2018).

Sedimentary pigments (0.2–0.5 g dry sample) were extracted using organic solvents of acetone: methanol: water (80:15:5, v/v/v) by leaving in a-20°C freezer for 24 h. Extracts were filtered with a 0.22 μm-pore PTFE filter, concentrated to dry with N₂ gas, and re-dissolved with organic solvents of acetone: ion-pairing reagent: methanol (70:25:5, v/v/v; Chen et al., 2001). Note that sedimentary pigments were measured at 0.5-cm resolution for upper 10 cm and 1-cm resolution for the lower. Pigments were measured via an Agilent 1200 series high performance liquid chromatography unit (HPLC, Agilent, United States) with quaternary pump, autosampler, 120EC-C18 Hypersil column (3.0 × 150 mm; 2.7 μm particle size), and photo-diode array detector (PAD). Pigments were identified and quantified based on the retention time and absorption spectra and authentic standards (DHI, Denmark; Leavitt and Hodgson, 2002; McGowan, 2013). The identified pigments included all algae (Chl α, pheophytin α), chlorophytes (Chl b, pheophytin b, lutein), siliceous algae (diatoxanthin), and cyanobacteria (canthaxanthin). The concentration of all pigments was expressed in nmol g TOC⁻¹. TOC, TN and TP were previously published in Chai et al. (2018). δ¹³C and δ¹⁵N were measured at 0.5-cm resolution using a Thermo Scientific MAT253 Delta V mass spectrometer coupled with a Flash elemental analyzer (Thermo Fisher, United States). The δ¹³C and δ¹⁵N is presented in δ notation as the per mil deviation (‰) from standards relative to Vienna PeeDee Belemnite (VPDB) and atmospheric N₂, respectively. Replicate analyses indicate that the analytical precision was better than 0.1‰ for the δ¹³C and 0.3‰ for the δ¹⁵N.

Principal components analysis (PCA) was used to summarise the major ecological trends of primary producer communities inferred from sedimentary pigment dataset. Redundancy analysis (RDA) was conducted to reveal the correlations between primary producer communities and multiple environmental factors. The response variables are sedimentary pigment concentration of the total algal community and specific algal group including chlorophytes, siliceous algae, and cyanobacteria. The explanatory variables include nutrient factors (TOC, TN, C/N, and TP), climate factors (temperature and precipitation), and atmospheric N deposition factors (δ¹⁵N and N₂O emission). Note that N₂O emission in the study region was obtained from grid monitoring data from EDGAR-Emissions Database for Global Atmospheric Research¹. In addition, variation partitioning analyses was further performed to investigate the synergistic effect of two main forcing factors on algal community dynamics. The climate changes (temperature) and atmospheric N deposition (δ¹⁵N) were selected as driving variables, while sedimentary pigment concentration of the total algal community and specific taxonomic groups as the response variables. All the variables were log-transformed prior to analysis, performed by the program Canoco 5 (Šmilauer and Lepš, 2014).

We further implemented a process-based lake ecosystem modeling approach for CQ Lake for a mechanistic understanding of the long-term ecological dynamics. Here we used the GOTM-WET model (Schnedler-Meyer et al., 2022), which included both General Ocean Turbulence Model (GOTM) and Water Ecosystem Tool (WET) for lake hydrodynamic and ecological modeling, respectively. GOTM is originally a 1D hydrodynamic model for ocean and its lake branch includes the hypsographic structure of the lake water and sediment, so that the model could simulate water temperature, ice cover, turbulence and mixing along the lake vertical dimension (Umlauf and Lemmin, 2005). WET is developed from the PCLake model (Janse et al., 2010), and can simulate the mass-balanced closed nutrient cycling and ecological processes with a customizable food web structure. Here we focused on the phytoplankton community and applied the ‘default’ setting that which includes three algae groups (diatoms, green algae, cyanobacteria) that are the dominant groups in CQ Lake. The model can simulate the lake dynamics at hourly basis. Considering that we primarily investigated the long-term lake ecosystem dynamic, we aggregated the model outputs from hourly to annually basis to better compare with the sedimentary data which are in annual time-resolution.

We collected the meteorological data (including wind speed, wind direction, air pressure, air temperature, relative humidity, cloud cover, solar radiation and precipitation) at hourly basis from the ERA5 global meteorological reanalysis data (from European Centre for Medium-Range Weather Forecasts, ECMWF) at the grid cell of CQ Lake (Supplementary Figure S2). In addition, we collected daily air temperature at the national meteorology station (CMA) ‘Shangrila’ (3,277 m a.s.l., ~20 km northeast of CQ Lake) from 1957 to 2020 AD, which is elevation-corrected with a reduction of ~5.2°C for CQ Lake. The Shangrila data is used for bias-correction of the ERA5 data (Supplementary Figure S3). Furthermore, the meteorological data without anthropogenic climate change were provided by the Inter-Sectoral Impact Model Inter-comparison Project (ISI-MIP) database as the ‘piControl’ scenario (Supplementary Figure S4), which were projected by four Global Climate Models (GCMs; Warszawski et al., 2014). The ISIMIP data were downscaled from daily to hourly basis following the method by Shatwell et al. (2019), except for wind speed and direction which were provided along the Eastward and Northward direction at sub-daily resolution (Frieler et al., 2017).

The annual discharge to the lake was estimated using the SCS-CN model (Kong et al., 2015) driven by annual precipitation data from ERA5. For nutrient loading from the drainage area (TN and TP), an export coefficient model (Johnes, 1996) was used considering vegetation (subalpine conifers and shrubs) as the land use type. We further incorporated annual atmospheric N deposition as part of the TN loading, based on the field observations in 1999 AD (2.22 kg hm⁻² yr.⁻¹; Zhang et al., 2013) and 2012 AD (20.0 kg hm⁻² yr.⁻¹; Zhu et al., 2015). The long-term data were estimated by the linear interpolation assuming a zero deposition in 1947 AD from the sediment N isotope peak (2.79‰; Supplementary Figure S5).

The hydrodynamic model (GOTM) was calibrated against the daily water temperature data from a neighboring Moon Lake at 1 m depth during 2013 (Supplementary Figure S6). The calibration was performed using the auto-calibration program PARSAC with several successful applications with GOTM (Ayala et al., 2020; Kong et al., 2022). For the ecological model WET, the initial parameter values were based on the default values calibrated to 40+ lakes (Janse et al., 2010), serving as a starting point for calibration. The initial values for state variables of nutrient (TN and TP) in the sediments were determined based on the sediment records at the corresponding layer.

We designed a model scenario analysis including both the TN loading and meteorological conditions as the two-factorial experimental simulations. For TN loading, there are two levels (N −/+) representing the scenarios of constantly low TN loading after 1980s AD, and the factual TN loading with an increase after 1980s AD and much rapidly after 2000s AD (Supplementary Figures S4, S5). For meteorological variables, two scenarios including the ‘piControl’ (PIC, or C−) and the ‘historical’ (HIST, or C+) were used representing the climate without and with anthropogenic impacts, respectively. In summary, the scenario analyses including four model scenarios (N −/+ and C −/+), while for C −/C+ there are four GCMs, resulting in 16 model runs in total.

Among the model runs, the scenario of TN+ and C+ represents the factual condition, which was considered as the ‘hindcast’ simulation. Due to lack of recent field observations for biogeochemical and ecological variables, we performed a model calibration during the hindcasting, by using the sediment records including the geochemical (TN and TP) and the pigment records of the three algal groups. A trial-and-error method was utilized to fit the model outputs from 1900 to 2020 AD to the sediment records. All the parameter values that were optimized during the calibration were documented in the Supplementary Table S1.

Next, we tested the hypotheses that the rapid increasing atmospheric N deposition and the climate change since the 1980s together dictated the ecological dynamics in CQ Lake, particularly in the succession of algal community. We evaluated if there are additive and/or synergistic effects of TN loading and climate change on the ecological dynamics. Relevant processes of the causal effects between the TN loading and climate change and ecological dynamics were provided by the process-based modelling with GOTM-WET. All statistical analyses were performed in R (R Core Team, 2021). Long-term trends and change rate in time series of model outputs were calculated based on the Mann-Kendall test and Sen’s slope by the R package ‘trend’ (Pohlert, 2020).

The chronology of CQ1 and CQ2 cores were previously reported in Chai et al. (2018) and Zhang et al. (2022), covering the past two centuries from 1800 to 2014 AD (Supplementary Figure S1). The sedimentary pigment records exhibited a two-stage pattern with relatively stable and low value before ~1980 AD and thereafter started distinctly increasing to the highest value in the top, expect for the cyanobacteria pigment (Figures 2, 3). In the interval of 1800–1980 AD, the pigments of siliceous algae (diatoxanthin) and cyanobacteria (canthaxanthin) were present in quite low abundances, while the pigments of total algae (Chl a and pheophytin a) and chlorophytes (Chl b and lutein) were relative moderate. Note that the concentrations of pheophytin a and pheophytin b are quite high and show a high value between 1800 and 1850 AD and shift to a low relative value afterward. In the interval of 1980–2014 AD, the total algae pigments began to increase rapidly until 2000 AD with a low value before. The chlorophytes pigments show a gradually increasing trend and an abrupt increase from the siliceous algae pigments, whereas the canthaxanthin pigment remains low in the whole profile.

For the pigment data, PCA axis 1 (PCA1) explained 71.9% of the total variance in the pigment assemblages and was positively correlated with all primary producers (Chl a and pheophytin a) and chlorophylls and their derivatives from chlorophytes (Chl b and pheophytin b; Figure 4A). PCA axis 2 (PCA2) explained a further 14.3% of the total variation and was positively correlated with cyanobacteria (canthaxanthin) and negatively correlated with chlorophytes (lutein) and siliceous algae (diatoxanthin). RDA results revealed that variations in pigment composition were significantly correlated with temperature and atmospheric N deposition (δ¹⁵N and N₂O). In contrast, other factors such as TN, TOC, TP, C/N, and precipitation show weak correlation with pigment composition changes (Figure 4B). The results of variation partitioning analyses reveal that climate warming and N deposition can independently explain 9.9 and 46.8% of the explained variables, respectively. In contrast, synergistic effects of climate warming and N deposition contribute anther 43.3% for the variability in algal community (Figure 4C). In contrast, climate warming, N deposition, and synergistic effects can explain 5.9, 27.7, and 25.6% of the total variables, respectively, with another 40.9% unexplained variables which may be due to the unconsidered additive or synergistic effects of multi-factors and the uncertainty of analytical method.

The ecosystem model GOTM-WET well captures the daily surface water temperature during 2013 (RMSE = 1.31°C, R = 0.92), with the seasonal variations of summer peaks (~15°C) and winter minimum (0°C; Supplementary Figure S6). The model also predicts an ice coverage from November to March, and the maximum ice thickness ranged 0.1–0.2 m, which are all in concert with field observations (Supplementary Figure S6). The model predicts that the lake stratifies from May to October with a mean thermocline depth of 3.11 m and a maximum of 9.49 m. The calibration during 2013 proves that the physic model (GOTM) is valid for the lake, so that the model could accurately predict the lake thermal dynamics for the years beyond 2013. The model hindcast from 1900 to 2020 AD provides the predictions of the lake physical dynamics driven by the factual climate and nutrient loading (Figure 5). The results suggest that the physic condition including water temperature, light in the epilimnion, and maximum ice thickness were rather stable before the 1970s. Light intensity started to increase from late 1970s AD, and temperature and ice started to change rapidly after the 2000s. From 2000 to 2020 AD, annual average water temperature increased at a rate of 0.017°C yr.⁻¹, while maximum ice thickness decreased at −0.07 cm yr.⁻¹. Annual average light in the epilimnion decreased at a rate of 0.1 W m⁻² yr.⁻¹.

The ecological model (WET) is calibrated during the hindcast against the long-term sediment records, in particular to the phytoplankton community indicated by the pigment data (Figure 6). The results show that the model nicely captures the long-term patterns of Chl a concentration and the biomass of two major algal groups (green and diatom). In specific, the rapid increase of Chl a and green algae since the 2000s are both represented, while the increase in diatom biomass since the 1980s are also well simulated. Though the absolute Chl a and biomass are difficult to quantify by the pigment data, the consistent patterns between pigment and model imply that the model has grasped the key succession features of the phytoplankton community over the historical period since 1900 AD. The model predicts the rather stable TN concentration in the water column before 2000 AD, and a rapid increase TN afterwards (Supplementary Figure S5), associated with the rapid increase TN atmospheric deposition. Taken together, the model calibration during 2013 AD and the hindcast since 1900 AD as described above provide evidence for reliable outcomes from the model scenario and mechanistic analyses below.

The model scenario analyses reveal that the increasing TN enrichment via atmospheric deposition and lake warming caused by climate change are the two main factors promoting the algal proliferation and the community succession, not only individually but also synergistically (Figure 6). Model results suggest that without the two factors (scenario ‘C-/N-’), algal proliferation would not occur, while the community would maintain a persistence and constant biomass of green algae and the secondary dominating diatom. With the TN enrichment alone (scenario ‘C-/N+’), Chl a is enhanced by 51%, green algae would flourish from the 2000s (+43%), while diatom also increases (+10%) but starts earlier around 1980s. This pattern is highly consistent with the sediment pigment data (Figure 6). This result implies that the increasing atmospheric N deposition acts as the key role in driving the phytoplankton succession, featured by the increasing biomass of green algae and diatom. Note that cyanobacteria abundance remains low in both sediment records and model predictions. On the other hand, under the lake warming alone (scenario ‘C+/N-’), however, the change is less visible compared to ‘C-/N+’, but still effective in the way that Chl a (+5%), green algae (+3%) and diatom (+2%) are slightly promoted in biomass. This result suggests that lake warming acts as a secondary role in driving the algae proliferation. Finally, when the two factors are simultaneously effective (scenario ‘C+/N+’), the phytoplankton community depicts the factual pattern as the historical succession, and the increases in Chl a (+77%), green algae (+67%), and diatom biomass (+19%) are even higher than the additive effects of each factor separately, indicating the existence of synergistic effects from both factors.

CQ Lake is located in the southeast of Tibetan Plateau, which belongs to one of the most sensitive areas to climate changes with twice as much warming as the global average (Immerzeel et al., 2020). Over the past century, meteorological data for CQ Lake from ECMWF reveals that this study region has experienced large-magnitude warming with a rapid increase of 1.5°C since 1980 AD (Figure 3A). Regional annual precipitation shows a gradual decrease trend between 1860 and 1990 AD, followed by an increase afterwards (Figure 3B), which can be supported by reconstructed regional precipitation and relative humidity records from tree ring and speleothem (An et al., 2013; Tan et al., 2016; Xu et al., 2019). Although the catchments of remote subalpine lakes in the Tibetan Plateau are rarely affected by anthropogenic activities directly, reactive N deposition rate in the mountains is high with a value of 20–40 kg N ha⁻¹ yr.⁻¹ (Liu et al., 2011; Zhang et al., 2012), which is mainly due to increased use of synthetic nitrogen fertilizer and fossil fuel combustion after 1980s (Liu et al., 2013). In the CQ lake, δ¹⁵N data display a distinct decline towards to 0‰ since 1980 AD (Figure 3C), which can be interpreted as a signal of enhanced atmospheric N deposition (Holtgrieve et al., 2011). Although δ¹⁵N in the sediments is affected by N inputs from different sources such as catchment sources, atmospheric deposition and in-lake cycling (Catalan et al., 2013), catchment vegetation and soil should contribute enrichment of δ¹⁵N in lake sediments (Bartrons et al., 2010). The effect of bacterial process on δ¹⁵N cannot explain the long-term significant depletion trend since 1980 AD (Talbot, 2001). Furthermore, monitoring data of global greenhouse gas emissions demonstrated a remarkable increase of N deposition at our study area inferred from N₂O concentration (Figure 3C). The consistent variations can also be observed in nearby Moon lakes with obvious depletion of δ¹⁵N (Hu et al., 2014) and in Himalayan ice cores with increased nitrate concentration (Thompson et al., 2000). Therefore, all lines of these evidence support that our study area has experienced unprecedented climate warming and high atmospheric N deposition over the past several decades.

Previous studies show TOC and TN can indicate the variation in primary productivity within the lake and catchment basin (Meyers and Ishiwatari, 1993; Chai et al., 2018). At 1800–1980 AD, the consistent decreasing TOC and TN indicates a long-term decline in productivity in the CQ Lake and its catchment (Figures 3E,F; Chai et al., 2018), which may be related to regional precipitation changes (Meyers and Ishiwatari, 1993; Zhang et al., 2019). Decreased precipitation may limit vegetation and soil development and further diminish the input of catchment organic matter into the lake (Chen et al., 2018), resulting in the gradual decline of TOC and TN. Similar trend can be observed in the total phosphorus of lake sediments (Figure 3D). After 1980 AD, the slight increase in TOC and TN indicates the increased inputs of allochthonous organic matter and/or enhanced primary productivity in the lake (Figures 3E,F; Chai et al., 2018). Nevertheless, the decoupling between total phosphorus and regional precipitation points to the main contribution from enhanced primary productivity rather than allochthonous organic matter inputs at this interval. In addition, evidence of enhanced primary productivity can be observed from the proxies of C/N ratio and δ¹³C (Figures 3G,H), because enhanced primary productivity within the lake can lead to decreased C/N ratio and depleted δ¹³C.

Cluster analysis of pigment assemblages in the CQ Lake reveals a remarkable shift in algal communities in ~1980 AD and ~ 2000 AD (Figures 3I–L, 4A). The concentration of most pigments showed a relatively stable low value before 1980 AD, and then siliceous algae started to increase to maximum value toward to present, and total algae and chlorophytes remained relatively stable followed by rapid proliferation until 2000 AD. In contrast, the concentration of cyanobacteria still kept stable and low value over the whole period (Figures 3I–L). Serving as dominant producers in aquatic ecosystems, algae can reflect variations of primary productivity and nutrient level (Reynolds, 2006). The remarkable increase in pigment indicates accelerating algae growth and enhanced primary productivity since 1980 AD. Although sedimentary pigment can be influenced by degradation effects downcore, the relatively stable ratio of Chl a and its degradation product Pheophytin a indicate stable preservation conditions in the CQ Lake (Figure 3M), similar to many other deep lakes (McGowan, 2013; Chen et al., 2018). The unambiguous increase in primary productivity and ecological changes can be widely observed in neighboring subalpine lakes from pigment and diatom proxies (Hu et al., 2014, 2017; Kong, L. et al., 2017; Wang et al., 2020). For example, the pigment PCA1 in nearby Shade Co shows consistent increase in primary productivity in the last few decades (Figure 3N), and similar trend can also be recorded in diatom composition (Figure 3O). Therefore, our multi-proxy paleolimnological records indicates that CQ Lake experienced remarkable environmental and ecological variations including intensified primary productivity and shift in algal community composition especially after 1980 AD.

Nutrient supply has been widely recognized as dominate drivers of lake ecosystem dynamics (McGowan, 2013; Yang et al., 2017). Our RDA results suggest that atmospheric N deposition (δ¹⁵N and N₂O emission) was strongly positively associated with the variations in primary productivity and algae abundance, while sediment phosphorus and TN show a slight negative correlation (Figure 4B). As a small subalpine lake without direct human disturbance, nutrient supply of CQ Lake mainly derives from catchment natural inputs of nitrogen and phosphorus, atmospheric N deposition, and altered in-lake nutrient cycling such as phosphorus sediment release (Hu et al., 2014). On the other land, the long-term decline in sediment phosphorus concentration plays a limited effect for algae proliferation especially since 1980 AD (Figure 3), although sediment phosphorus cannot be completely equivalent to soluble phosphorus in lake, which may be affected by sediment release and atmospheric deposition (Catalan, 2000). Modern monitoring shows atmospheric phosphorus deposition is very low compared with regional atmospheric N deposition (Wang et al., 2017). Although phosphorus can be release from sediments and catchment weathering, the amounts is generally neglectable comparing to catchment inputs in small oligotrophic lakes such as CQ Lake (Catalan, 2000).

Assuming that sediment phosphorus is dominated by catchment inputs, the nitrogen content form catchment should show a similarly reduced trend. Whereas the slightly increase in sediment TN after 1980 AD in the CQ Lake mainly reflects the contributions of other sources such as atmospheric deposition and in-lake cycling (Figure 3; Thompson et al., 2000). Previous study suggests that increased N deposition is indeed an important factor for changes in nutrient sources (Hobbs et al., 2011; Hu et al., 2014). Our model results supported a long-term increase of TN loading after 1980 AD (Figure 3F) when accounting for additional atmospheric deposition inputs even under the condition of reducing catchment inputs (Supplementary Figure S5). In addition, atmospheric N deposition can stimulate increase in primary productivity especially in subalpine lakes. Further quantitative investigation reveals that N deposition with more than 10 kg ha⁻¹ yr.⁻¹ would cause significant ecological changes for the remote lakes (Bergstrom and Jansson, 2006). In the southeast of Tibetan Plateau, the deposition rate is 20–40 kg N ha⁻¹ yr.⁻¹, which is close to or above the estimated critical load of N deposition for ecological changes (Liu et al., 2011; Zhang et al., 2012). At present, CQ Lake has relatively high phosphorus and low nitrogen concentration with TN:TP of ~8 and belongs to a N-limited lake. The N-limited condition may be more significantly considering the gradually declining sediment phosphorus and increasing nitrogen concentration over the past several decades (Figure 3).

Under N-limited conditions, remote lakes are quite susceptible to the enrichment effect of atmospheric N deposition, especially when deposition rates exceed the critical load (Baron et al., 2011). After1980 AD, increased active N deposition largely promotes primary productivity in the CQ Lake charactered by increase in the abundance of total algae, diatom and chlorophytes (Figure 3). Note that stable concentration of canthaxanthin pigment at this interval shows that the lake nutrient level remains too low to support the growth of cyanobacteria. Furthermore, model scenario analysis suggests nutrient loading is indeed critical factor to the variations of primary productivity and algal abundance in the CQ Lake (Jeppesen et al., 2005; Kong X. et al., 2017). Enhanced nutrient supply, caused by the rapid increase of atmospheric N deposition, stimulates algae proliferation characterized by remarkable increase in the biomass of Chl a, green algae, and diatom despite the changes in climate conditions (Figure 6). Likewise, enhanced primary productivity over the past several decades can also be observed in neighboring other lakes such as Jiren Lake (29°43′ N, 100°48′ E), Tiancai Lake (26°38′ N, 99°43′ E) and Wuxu Lake (29°09′ N, 101°24′ E; Wischnewski et al., 2011; Hu et al., 2017). Overall, our results further approve that the influence of atmospheric N deposition on alpine lakes is magnificent in the Tibetan region.

Our RDA results suggest that climate warming also play a positive role in algae abundance and community succession in the CQ Lake (Figure 4). Climate warming can stimulate changes in primary producer communities by direct physics process such as heating surface lake water, shortening the duration of ice-cover, and altering thermal stratification (Catalan et al., 2013; Woolway et al., 2020). Generally, lake warming and enhanced lake stratification due to climate changes provide more stable and suitable conditions for the growth of primary producers (Hébert et al., 2021), and prolonged growing season due to shortened ice-cover periods also provides more time for algae proliferation (Rühland et al., 2015). However, our model results reveal a limited sole contribution of climate warming to the increase in lake primary productivity and algal proliferation over the past several decades in the CQ Lake (Figure 6). The low contribution may be due to relative low warming extent by ~0.6°C increase in water temperature since 1980 AD (Figure 5A). This limited variable amplitude is insufficient to trigger detectable changes in the ecological process (Hébert et al., 2021; Kong et al., 2021). Nevertheless, our model shows that the effects of climate warming on lake primary productivity and algae proliferation are more pronounced under increasing atmospheric N deposition (N+; Figure 6), suggesting a synergistic effect of climate warming and N deposition on lake ecological changes. Meanwhile, variation partitioning analyses was performed to investigate the synergistic effect of two main forcing factors on algal community succession. The results suggest that climate warming and N deposition can independently explain 9.9% and 46.8% of the explain variables respectively, while their synergistic effects contribute anther 43.3% for the algae growth and community succession (Figure 4C). Such agreement in analytic results from two independent approaches (statistical analysis and model simulation) points to the reliability of the findings that the major role of atmospheric N deposition and the synergy with lake warming driving the long-term ecological dynamics in the CQ Lake.

Furthermore, the process-base modeling further unravels the potential mechanisms that may formulate such synergy. On the one hand, lake warming may increase the nutrient availability due to higher flux of temperature-driven mineralization, which may exacerbate the impacts of excessive nutrient after atmospheric deposition especially after 2000s AD. Similarly, in situ experiments in Rocky Mountain lakes revealed that small planktonic species respond to lake thermal structure changes only with adequate nutrient availability (Saros et al., 2012), suggesting that the effect of lake warming depends on the nutrient availability that cycling between dissolved and particulate forms. On the other hand, the declining maximum ice thickness in the lake since 1980 AD due to warming may increase the light availability under ice (Figures 5B,C), which may act as another indirect effect of lake warming and further promote the synergistic effects with excessive nutrient. Previous studies have implied that global warming leads to mild winter and increased light intensity related with thinner ice coverage, which may promote algae proliferation via more intensive photosynthesis and higher efficiency of nutrient utilization (Kong et al., 2021). Furthermore, increasing light intensity in water may also promote algae proliferation via more intensive photosynthesis and expansion of benthic habitat. Thus, winter is not a dormant period and may even dictate ecological processes in other seasons (Hébert et al., 2021; Hrycik et al., 2021). In addition, our sediment pigment also indicates an earlier diatom growth than other algae since 1980 AD (Figures 3I–L), and this can also be observed in the model result (Figure 6C), which is in coincidence with a rapid increase in light intensity (Figure 5B). This is because diatom is highly sensitive to light changes in lake. The ‘light switch’ after lake stratification has been conventionally considered as the trigger of spring diatom blooms (Sommer et al., 2012). We further hypothesize that the increasing in light after 1980 AD was because of thinner ice coverage driven by lake warming, which further enhance the synergy with nutrients. However, such hypothesis needs to be tested in future studies with more reliable field data on ice and light availability from the lake.