Synergistic effects of climate warming and atmospheric nutrient deposition on the alpine lake ecosystem in the south-eastern Tibetan Plateau during the Anthropocene

1. Introduction

The Tibetan Plateau, as the “Third Pole” of the Earth (Qiu, 2008), is one of the most sensitive areas to climate changes at a rate twice warming than the global average from meteorological observations (Immerzeel et al., 2020). Lakes are sentinel of global changes and sensitive to unprecedented climate warming and pervasive human activities (Adrian et al., 2009). Evidence from paleolimnological investigations suggests that aquatic ecosystems of montane lakes have correspondingly experienced rapid and pronounce changes in recent decades on the plateau (Hu et al., 2014, 2017). Climate warming can stimulate changes in primary producer communities through directly heating surface lake water, shortening the duration, thinning the ice-cover, and altering thermal stratification, or indirectly affects ecological changes by interaction with nutrient recycling (Catalan et al., 2013; Woolway et al., 2020). Recent studies indicate algal communities and diatom compositions display a trend towards long-term warming over the past century (Hu et al., 2017; Kong, L. et al., 2017). Specifically, reduced ice cover and enhanced thermal stratification caused by climate warming induce the diatom shifts from benthic Fragilaria taxa to planktonic Cyclotella with a preference of warm condition (Smol et al., 2005; Hu et al., 2014; Rühland et al., 2015).

Anthropogenic impacts in populated area are migrating to these remote areas via transporting particulate nutrients by atmospheric deposition (Wolfe et al., 2003). Tibetan Plateau region is now a hotpot of atmospheric nitrogen (N) deposition, transporting from lower latitude southeast Asia by monsoon circulation, with the N deposition rate of up to 20–40 kg ha⁻¹ yr.⁻¹ (Chen et al., 2004; Zhang et al., 2012). Furthermore, nitrate concentration recorded in Himalayan ice cores revealed a long-term increase in N pollution since 1850 AD (Thompson et al., 2000). Previous study suggests that increased N deposition is indeed an important nutrient source, promoting primary productivity and driving ecological changes when N deposition rate exceeds ~10 kg ha⁻¹ yr.⁻¹ at remote montane lakes (Bergstrom and Jansson, 2006; Hu et al., 2014). Such a high N deposition rate in the southeast of Tibetan Plateau should be close to or above the estimated critical load for ecological changes (Liu et al., 2011; Zhang et al., 2012). Furthermore, lake primary productivity enhancement and phytoplankton composition alternation driven by N deposition are also reported especially for N-limited alpine lakes in the Tibetan Plateau region (Hu et al., 2014).

Ecological dynamics of aquatic system are not only driven by a single stressor such as climate change, human activity, or nutrient supply, but more often by multiple stressors (Birk et al., 2020), which have been increasingly recognized as the joint drivers of lake ecological changes (Ormerod et al., 2010). A recent study from 33 freshwater mesocosm experiments indicates that more than 60% of the cases are subject to additive effects and interactive (synergistic or reversal) effects by the paired-stressor combinations (Birk et al., 2020). Furthermore, a modern survey based on more than 1,000 lakes in U.S. suggests that the synergistic effects of climate changes and nutrients supply are indeed widespread in determining lake primary productivity and algal community shifts, especially for eutrophic lakes (Rigosi et al., 2014). In addition, in situ experiments in Rocky Mountain lakes also reveal small planktonic species respond to lake thermal structure changes only with adequate nutrient availability, which also supports the interactive effects of lake warming and nutrient availability for mesotrophic-oligotrophic alpine lakes (Saros et al., 2012). However, whether warming and atmospheric deposition would impose additive or even synergistic effects on the lake ecosystem in the Tibetan Plateau region are rarely investigated and remains unknown. Therefore, in-depth research on the synergistic effects of multiple stressors is warranted for improved understanding the response mechanisms of lake ecosystem on the Tibetan Plateau.

The long-term perspective is crucial for understanding the nature and timing of the lake ecological evolution, and its feedback mechanism responding to external multiple stressors (Smol, 2010). Paleolimnological approaches are widely used to reconstruct the long-term dynamics of environmental changes and ecological evolution based on multi-proxy data over decadal to centennial timescales (Zhang et al., 2019; Fordham et al., 2020), exemplified by numerous studies in European, North American, Arctic, and Tibetan Plateau region (Smol et al., 2005; Rühland et al., 2015). For example, fossil pigment proxy can reflect the pelagic algal community succession of the lake, and the tempting outcomes from sediment pigment shed lights on exploring the phytoplankton community exemplified from our study and others in more eutrophic lakes (Leavitt et al., 1997; Leavitt and Hodgson, 2002; McGowan, 2013). However, the key drivers and underlying mechanisms of such dynamics remains difficult to quantify based on the sediment records alone. On the other hand, process-based ecosystem modelling offers the opportunity to hindcast the lake evolution history and to identify the key drivers (climate and human activities). Yet such approach is constantly restricted by lack of external forcing (climate, nutrient loading etc.,) and field data, particularly in those understudied lakes without a routine monitoring program (Kong, X. et al., 2017). Combining the modelling and paleolimnological approach is therefore promising to improve our understanding of the long-term lake evolution, which calls for more investigations including the Tibetan Plateau region (Seddon et al., 2014).

In this study, a representative subalpine lake (Cuoqia Lake) in the southeast of Tibetan Plateau was chosen as a monitor of regional environmental and ecological changes. We use high-resolution paleolimnological records to (1) investigate the regional environmental changes and lake ecological evolution (primary productivity and algal biomass) over the past two centuries based on ²¹⁰Pb/¹³⁷Cs dated multiproxy data including geochemistry, stable isotopes and sedimentary pigments, (2) disentangle its response mechanisms to climate warming and atmospheric nutrient deposition using ecosystem modelling, and (3) explore the synergistic effects of multiple external stressors in regulating primary productivity and algae abundance changes through statistical analyses and comparative scenario modelling.

2. Study site

Cuoqia Lake (CQ Lake, 27°24′ N, 99°46′ E, 3960 m a.s.l.) is located in the hinderland of Hengduan Mountains in the southeast margin of Tibetan Plateau, ~20 km southwest of Shangri-la country, Yunnan province (Figure 1A). It is a small subalpine ice-scour lake with the maximum water depth and average water depth of ~26.4 m and ~ 13.2 m, respectively. CQ Lake have a surface area of ~0.07 km² and a catchment area of ~0.44 km² with a small catchment and lake area ratio of 6.3 (Figures 1B,C). It is a semi-closed lake without visible surface inflows and outflows, and is mainly fed by seasonal rainfall and snowmelt water. The lake is surrounded by forest vegetation on three sides and steep rock walls on one side. There are no visible aquatic plants in lake water with a pH of ~6.7 and dissolved organic carbon (DOC) of 5.15 mg/l in May 2014. The watercolor is close to sauce-brown with a Secchi depth of 3.7 m. Field survey indicates that it is an oligo-trophic lake with total phosphorus (TP) and total nitrogen (TN) of ~10 μg L⁻¹ and ~80 μg L⁻¹, and is N-limited with TN:TP of ~8 (Hu et al., 2017; Li et al., 2021).

Figure 1. Location and setting. (A) Map showing the location of the study site Cuoqia Lake (CQ Lake, red dots) and nearby Moon Lake and Shade Co (gray dots) in the southeast Tibetan Plateau. (B) topography and catchment (circled by yellow line) of CQ Lake, base map from Google Earth [modified from Li et al. (2021)]. (C) Bathymetry of the CQ Lake and locations of the sediment cores (red asterisks). (D) Mean monthly temperature and precipitation at the Shangri-la Meteorological Station, Yunnan Province, China (data retrieved from the Meteorological Administration of China, http://data.cma.cn/).

Currently, the study site belongs to highland humid monsoon climate region, and has distinct rainy and dry seasons with wet in summer and dry in winter. It is mainly influenced by the interaction of Indian summer monsoon, the southern branch of the Westerly wind, and local Tibetan Plateau climate (Zhang et al., 2022). There is no weather station around the lake, and the nearest Shangri-la meteorological station (3276.7 m a.s.l.) shows that the mean annual precipitation (MAP) is ~651.1 mm with more than 81.3% of precipitation falling during monsoon season from April to September over the period of 1958–2020. The mean annual temperature (MAT) ~6.3°C with a maximum monthly mean of ~13.9°C in July and a minimum of ~ − 2.3°C in January, and the relative humidity is ~68% (Figure 1D). At present, the lake lies below the tree line and surrounded by forest vegetation on three sides and steep rock walls on one side. Primary vegetation around the lake appears to be undisturbed in the historic past, dominated by cold temperate conifer forest such as Abies and Picea, shrub such as azalea, Yunnan wormwood such as Kobresia yunnanensis and miscellaneous grass meadow (Xiao et al., 2014; Li et al., 2021).

3. Materials and methods

3.1. Sample collection and age control

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).

3.2. Proxy analysis

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.

3.3. Data analysis

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).

3.4. Lake ecosystem modelling configuration

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.

3.5. Model hindcast and scenario analyses

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).

4. Results

4.1. Sedimentary records

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.

Figure 2. Sediment pigment. The concentration (nmol g⁻¹ TOC) of main pigment including (A) Chl α, (B) pheophytin α, (C) Chl b, (D) lutein, (E) pheophytin b, (F) diatoxanthin, and (G) canthaxanthin in the sediments at the CQ Lake.

Figure 3. Historical environmental and ecological changes around the CQ Lake. (A) temperature and (B) precipitation changes, (C) sedimentary δ¹⁵N and N₂O emission in the study region, (D) the concentration of sediment TP (Chai et al., 2018), (E) %TOC (Chai et al., 2018), (F) %TN (Chai et al., 2018), (G) the ratio of C/N (Chai et al., 2018), (H) δ¹³C; The biomass of different algae (average concentration from sediment pigments) including (I) total algae (Chl α and pheophytin α), (J) chlorophytes (Chl b, lutein, and pheophytin b), (K) siliceous algae (diatoxanthin), (L) cyanobacteria (canthaxanthin); (M) The ratio of Chl α/Pheo α, (N) pigment-based PCA at the Shade Co, and (O) diatom-based DCA at the Lugu Lake.

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.

Figure 4. (A) Principal components analysis of pigments and (B) redundancy analysis between pigments and significant explanatory variables including temperature, precipitation, δ¹⁵N, N₂O emission, TOC, TN, C/N, and TP. (C) Variation partitioning analyses was further performed to investigate the synergistic effects of two main forcing factors (climate warming and atmospheric nutrient deposition) on algal community dynamics.

4.2. Model calibration and hindcast

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.⁻¹.

Figure 5. Model simulated annual averaged data for (A) epilimnetic water temperature volume-weighted across the epilimnion, (B) light intensity in the lake epilimnion weighted by the area of the water–sediment interface in each modeled layer, (C) maximum ice thickness from 1900 to 2020 under different scenarios. The scenarios are composed for a two-factorial approach. ‘C’ denotes the climate factor (‘C-’ is the ‘piControl’ climate scenario without anthropogenic change, and ‘C+’ is the ‘historical’ with climate change). ‘N’ denotes the nutrient loading factor (‘N-/N+’ represents the scenarios without and with the rapid increase of atmospheric N deposition after ca. 2000s). The shade area represents the modeling range derived from the four Global Climate Models (GCMs) from the ISIMIP project. (D-F) The corresponding boxplots of the variable above during 2000 to 2020 under various modeling scenarios.

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.

Figure 6. Model simulated annual averaged, volume-weighted data across the water column for (A) Chl a concentration (B) biomass of green algae, (C) biomass of diatom, from 1900 to 2020 under different scenarios. The scenarios are composed for a two-factorial approach. ‘C’ denotes the climate factor (‘C-’ is the ‘piControl’ climate scenario without anthropogenic change, and ‘C+’ is the ‘historical’ with climate change). ‘N’ denotes the nutrient loading factor (‘N-/N+’ represents the scenarios without and with the rapid increase of atmospheric N deposition after ca. 2000s). The shade area represents the modeling range derived from the four Global Climate Models (GCMs) from the ISIMIP project. (D–F) The corresponding boxplots of the variable above during 2000–2020 under various modeling scenarios.

4.3. Model scenario analyses

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.

5. Discussion

5.1. Ecological and environmental changes at the CQ Lake over the past two centuries

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.

5.2. Nutrient deposition driving ecological evolution at the CQ Lake

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.

5.3. Synergistic effects of climate warming and nutrient deposition

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.

6. Conclusion

Our high-resolution paleolimnological records reveal the environmental and ecological variations at the CQ Lake in the southeast Tibetan Plateau over the past two centuries based on ²¹⁰Pb/¹³⁷Cs dated multiproxy data, including geochemistry, stable isotopes and sedimentary pigments. The most remarkable ecological changes indicated by sedimentary pigments at the CQ Lake have occurred since 1980 AD, manifested by a substantial increase in primary productivity and algae biomass. Our RDA results suggest that atmospheric nutrient deposition and climate warming were positively associated with the variations in primary productivity and algae abundance, indicating the underlying driver-response relation between both factors.

Ecosystem modelling, established based on contemporary data and paleolimnological records, was applied to disentangle the driving mechanisms of lake ecosystem dynamics. Scenario analysis revealed that nutrient deposition is the critical factor in driving the increase in lake primary productivity and algal proliferation over the past several decades in the CQ Lake, whereas the contribution of climate warming is limited under low nutrient conditions. Nevertheless, the effect of climate warming on lake ecological changes exhibits more pronounced effects under increased nutrient supply, indicating a synergistic effect of climate warming and nutrient deposition, which can also be supported by the results of the variation partitioning analyses with 43.3% contribution for synergistic effect of climate warming and N deposition. The modeling approach provides the unique opportunity to interpret the potential mechanisms shaping the synergistic effects between lake warming and atmospheric nutrient deposition, which was conventionally difficult with data analyses alone. Our results highlight that multiple environmental stressors under global changes could not only additively, but also synergistically, affect lake ecosystems, which may possibly worsen the negative impact on remote, subalpine lakes as CQ Lake. Thus, more attention on such synergy is warranted to evaluate the impacts of global changes on aquatic ecosystems.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

CZ and XK conceptualization, data analysis, preparation of figures and tables, and writing. CZ, QL, and KZ developing methods. CZ, BX, LC, KZ, XY, and JS conducting the research and data interpretation. All authors contributed to the article and approved the submitted version.

Funding

This study is supported by the National Natural Science Foundation of China (42007401 and 42177062), the National Key Research and Development Program of China (2019YFA0607100), the Strategic Priority Research Program of Chinese Academy of Sciences (XDB40000000), and the Science and Technology Planning Project of NIGLAS (NIGLAS2022TJ02, NIGLAS2022GS02 and NIGLAS2022GS09).

Acknowledgments

We thank Yifan Chai, Lingyang Kong, Qian Wang for field and laboratory assistance.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1119840/full#supplementary-material

Footnotes

References