Using Wavelet Analysis to Examine Long-Term Variability of Phytoplankton Biomass in the Tropical, Saline Lake Alchichica, Mexico

: The phytoplankton biomass (chlorophyll-a, Chl-a) is directly related to the total production of lakes. Chl-a in temperate lakes oscillates on an annual scale. However, Chl-a oscillations in tropical lakes have hardly been documented, particularly over multiple years. Here, we described the periodicity of the Chl-a by performing a continuous wavelet analysis of 21 years (1998–2018), monthly Chl-a data from tropical, saline Lake Alchichica, Mexico. Parallel wavelet analyses were made on environmental time series (i.e., euphotic zone, mixed layer, temperature, dissolved oxygen concentration, dissolved inorganic nitrogen, soluble reactive phosphorus, soluble reactive silica). Throughout the time series, the wavelet transforms identiﬁed a regular and predictable annual cycle of the Chl-a associated with the warm-monomictic thermal-mixing pattern, the variability of the annual Chl-a cycle, and the presence of other cyclicities, 2-year and ~4–5 years, associated with external forcing agents (e.g., North Paciﬁc Oscillation). The water quality variables display a recurrent annual cycle. At the same time, the trophic variables (nutrient concentration) showed the same cyclicity as Chl-a (1-year, 2-year, and 4-year), suggesting the external forcing agents promote Chl-a augment through nutrient increase made available from stronger, deeper, mixing periods. Chl-a cycle showing large inter-annual variability, (c) identiﬁed the presence of a biennial cycle most likely controlled through SRSi availability and Si/N ratio associated to climate variability inﬂuencing the length and strength of the annual thermal stratiﬁcation, and (d) recognize other cycles with different periodicities associated with external forcing agents (e.g., ENSO-El Niño Southern Oscillation). The water quality variables were regular along the 21 years, displaying a recurrent annual pattern. Nutrient concentration, as trophic indicators, revealed the same cyclicity (1-year, 2-year, and ~4–5 years) as Chl-a, suggesting the external forcing agents lead to deeper, more complete mixing events increasing nutrient availability, which in turn stimulate larger Chl-a concentrations.


Introduction
Phytoplankton biomass, expressed as chlorophyll-a (Chl-a), is directly related to the total production of lakes [1] and can be used as a proxy for primary production [2]. Spatial and temporal changes in Chl-a are primarily controlled by light, temperature, and nutrients in the water column. These factors also determine the phytoplankton's size spectrum and species composition [3,4]. The annual phytoplankton succession was well documented in deep, thermally stratified temperate lakes (the Plankton Ecology Group model, e.g., [5,6]), with a unimodal phytoplankton maximum in the spring in oligotrophic lakes due to nutrient limitation.
Although many studies examined the phytoplankton variability within years, and more recently, a few long-term phytoplankton data are becoming available, few studies encompass datasets of several years or decades, primarily due to the lack of data [7]. The absence of long-term phytoplankton data and lack of continuity in records is common for tropical lakes [8,9]. According to Talling [10], the annual patterns of phytoplankton seasonality in the tropical belt are dominated by hydrological features (water inputsoutputs) or by hydrographic ones (water-column structure and circulation). Lewis [11] mentioned that in tropical lakes, the changes in phytoplankton abundance could be driven autogenically, as the organisms modify their environment (e.g., through nutrient uptake) or allogenically, associated with abiotic influences (e.g., water column mixing). Phytoplankton biomass in large tropical lakes varies at different time scales from seasons to centuries [12]. The seasonal scale in the tropics is driven by the alternation between the rainy (warm) and dry (cold) periods. In contrast, the global climate (e.g., El Niño Southern Oscillation, the Pacific Decadal Oscillation) drives the long-term scale. Despite the apparent "endless summer" [13] in tropical inland waters, there is commonly one event of deep watercolumn mixing, but also other different time-variability patterns (period/frequency and amplitude/range) such as the diel (24 h) cycle, the annual cycles with other periodicities (lunar, seiches, sunspots), and long-term and aperiodic or non-seasonal environmental changes [4,14].
Annual phytoplankton patterns differ across ecosystems, have large year-to-year variability [15][16][17], and are not strongly expressed in all aquatic ecosystems [18,19]. Whereas high variability of seasonal phytoplankton patterns is documented, no systematic analysis of annual cycles has been conducted to identify their characteristic periods of biomass variability and the recurrence strength at those periods [20].
Wavelet analysis performs a time-frequency analysis of the signal, which permits the estimation of the spectral characteristics of the signal as a function of time [21] and then the identification of different periodic components and their time evolution all along with the time series. Wavelet transforms are useful tools for determining the dominant scales of variation in time series and are used to examine periodicities in plankton communities [15,18,20]. However, this approach has been scarcely applied in tropical limnology.
In this work, our main goal was to reveal the variability in the annual periodicity of the phytoplankton Chl-a over 21 years and relate the variability with the environmental variables in the tropical saline Lake Alchichica, Puebla, Mexico. For this purpose, we addressed the following research questions: What are the cyclicity periods throughout 21 years (1998-2018) of (1) water quality and trophic status and (2) Phytoplankton Chl-a?
Our working hypothesis was that in this tropical, saline lake, the warm-monomictic thermal mixing pattern is the key factor determining the phenology of the pelagic Chl-a concentration. Specifically, we expected (1) to identify a predictable and recurrent annual cycle and (2) to unravel longer-term cyclicity. In order to answer the research questions and test the proposed hypothesis, our approach was to evaluate the temporal, long-term changes using data derived from monthly monitoring of the water column for 21 years  in the (a) water quality and trophic status and (b) phytoplankton Chl-a.

Study Area
Lake Alchichica is a volcanic maar in Puebla State (19 • 24.7 N, 97 • 24.0 W, 2300 m above the sea level), at the easternmost portion of the Mexican Plateau [22] (Figure 1).
The climate is temperate and semi-arid with an average temperature of 12.9 • C, rainfall < 500 mm year −1 , concentrated in the ("summer") rainy season, and evapotranspiration rate up to 1690 mm year −1 [3]. It is a deep lake (maximum depth 62 m, mean depth 40.9 m), circular in shape (maximum length 1902 m), with a surface area of 2.367 km 2 , a shoreline of 5.84 km, and a volume of 114,688,900 m 3 [23]. The lake is saline (total dissolved solids = 10.08 g L −1 , electrical conductivity = 13.53 mS cm −1 ) with chloride, bicarbonate, sodium, and magnesium as dominant ions and an alkaline pH (8.7-9.2) [22,24]. Alchichica is a warm monomictic lake that circulates during the cold, dry season, from the end of December or the beginning of January to the end of March or the beginning of April. The lake remains stratified throughout the warm rainy season, and the hypolimnion becomes anoxic for most of this period [22]. Lake Alchichica is oligotrophic with an overall nutrient concentration of <4 µmol L −1 at the mixed layer [25]. The climate is temperate and semi-arid with an average temperature of 12.9 °C, rainfall < 500 mm year −1 , concentrated in the ("summer") rainy season, and evapotranspiration rate up to 1690 mm year −1 [3]. It is a deep lake (maximum depth 62 m, mean depth 40.9 m), circular in shape (maximum length 1902 m), with a surface area of 2.367 km 2 , a shoreline of 5.84 km, and a volume of 114,688,900 m 3 [23]. The lake is saline (total dissolved solids = 10.08 g L −1 , electrical conductivity = 13.53 mS cm −1 ) with chloride, bicarbonate, sodium, and magnesium as dominant ions and an alkaline pH (8.7-9.2) [22,24]. Alchichica is a warm monomictic lake that circulates during the cold, dry season, from the end of December or the beginning of January to the end of March or the beginning of April. The lake remains stratified throughout the warm rainy season, and the hypolimnion becomes anoxic for most of this period [22]. Lake Alchichica is oligotrophic with an overall nutrient concentration of <4 µmol L −1 at the mixed layer [25].
Nineteen species compose the phytoplankton community in Lake Alchichica, with diatoms the most diverse group with ten species. The annual phytoplankton dynamics Nineteen species compose the phytoplankton community in Lake Alchichica, with diatoms the most diverse group with ten species. The annual phytoplankton dynamics display three characteristic features: (1) a winter diatom bloom during the mixing season, (2) a cyanobacterial bloom at the onset of the stratification period, and (3) a deep chlorophyll maximum (DCM) along the stratification period [26].
The phytoplankton biomass in Lake Alchichica is dominated (>70%) by large-sized individuals [3] [27] comprised most of the phytoplankton biomass. It dominated the winter diatom bloom and the DCM during the stratification period, except for a brief period at the onset of the stratification when Nodularia aff. spumigena bloomed [28].

Physical, Chemical, and Chl-a Characterisation
The long-term database herein analyzed comprises a monthly sampling program made over 21 years (1998-2018) for Chl-a, euphotic zone, temperature, and dissolved oxygen in the mixed layer, and over 20 years (1999-2018) for nutrients. All samplings were made at the same station, located in the central and deepest part of the lake, and at the same time (12-13 h).
Vertical photosynthetically active radiation (PAR) and natural fluorescence profiles were recorded with a Biospherical natural fluorescence profiler model PNF-300. The profiles consisted of discrete readings recorded every second as the sensor was lowered at about 1 m per 3-4 s, providing a vertical resolution of 25-30 cm. PAR profiles were used to identify the depth of the euphotic zone (Z EU = 0.1% surface PAR). Vertical water temperature and dissolved oxygen profiles were recorded with water quality monitoring probes (Hydrolab Datasonde models DS3, DS4, or DS5) with a vertical resolution of 1 m. Temperature and dissolved oxygen (DO) profiles were used to recognize the main stages of the lake's thermal regime (stratification and mixing) and to identify the depth of the mixed layer (Z MIX ).
Chlorophyll-a concentration (mg m −3 ) was calculated from the natural fluorescence flux, F f , and the incident irradiance. The natural fluorescence of Chl-a is defined as the total flux of light emitted by Chl-a in a suspension of phytoplankton per unit volume under ambient light. Keifer et al. [29] and Chamberlain et al. [30] showed that measurements of natural fluorescence flux (F f ) and scalar irradiance of PAR (E o (PAR)) provide accurate estimates of Chl-a. Chl The natural fluorescence records of the first five surface meters were omitted since the algorithms (from the profiler software) used to calculate the concentration of Chl-a from the natural fluorescence interpret that all the flow of red light directed upwards comes from Chl-a fluorescence. However, near the water's surface (the first five meters in transparent waters such as Alchichica), this interpretation is incorrect due to the dispersion of the red wavelength [29].
There are two important optical assumptions used in this: o a c (PAR) is the chlorophyllspecific absorption coefficient (absorption normalized to chlorophyll concentration), and ϕ f is the quantum yield of fluorescence. These values are treated as constants in the software for the PNF-300 with values of 0.04 m 2 mg −1 and 0.045 mE fluoresced per mE absorbed, respectively. The Chl-a concentration profiles were integrated per unit area (m 2 ) to calculate the total Chl-a concentration per m 2 .
Before wavelet analysis, we performed a seasonally-adjusted Mann-Kendall Trend Test to account for any seasonality in the phytoplanktonic chlorophyll-a concentration , water quality (1998-2018), and trophic (1999-2018) variables of Lake Alchichica. Analyses were run in R-4.0.3 [32].

Wavelet Analysis
Wavelet analysis used continuous wavelet transforms to determine the dominant scales of variability for Chl-a, physical and chemical time series from 1998 to 2018 as a function of period and time. Wavelet analysis emerged as a tool for characterizing periodicities in non-stationary time series, as it decomposes a time series in both the frequency and time domains [33]. The latter is why it has become an important tool to analyze time series since it performs a time-frequency analysis of a signal, allowing the estimation of the spectral characteristics of the signal as a function of time [21,34], and at the same time, identifies the different periodic components and their time evolution over time. The use of wavelet analysis in aquatic ecology has increased over the last few years, especially in marine ecology. Wavelet transforms are useful tools for determining the dominant scales of variation in time series and are increasingly used to examine the periodicities in plankton (e.g., [15,20]) and other aquatic communities. Wavelet analysis performs a time-scale decomposition of the signal by estimating its frequency characteristics as a function of time [20,35].
The continuous Morlet wavelet transform was used to normalize monthly time series data and obtain the wavelet power spectrum of the time series. The wavelet power spectrum identifies the periods that are the most important sources of variability across time. In addition, it is possible to obtain a global wavelet spectrum, which is similar to Fourier spectra, that identifies the variance associated with each period for a given time series [33]. Due to the nature of the data, the original time series is not equally spaced and was made regular using spline interpolation. All the analyses for spline interpolations and wavelet transforms were performed using Matlab.

Chlorophyll-a
We observed substantial variability in the Chl-a concentration in Lake Alchichica  The Chl-a concentration fluctuated interannually ( Figure 3). For 48% of the period, the Chl-a concentration was above the annual average, while 52% was below the annual average. The above and below the average Chl-a concentration regularly alternated (three above, two below, three above, three below, two above, four below, two above, two below). Annually, Chl-a concentration displayed a regular dynamic with two peaks ( Figure  3). The first peak was associated with the winter (January to March) mixing period when diatoms bloom throughout the water column. The second peak was related to developing a DCM at the metalimnion during the stratification reaching maximum concentrations in September-October. The Chl-a concentration fluctuated interannually ( Figure 3). For 48% of the period, the Chl-a concentration was above the annual average, while 52% was below the annual average. The above and below the average Chl-a concentration regularly alternated (three above, two below, three above, three below, two above, four below, two above, two below). Annually, Chl-a concentration displayed a regular dynamic with two peaks (Figure 3). The first peak was associated with the winter (January to March) mixing period when diatoms bloom throughout the water column. The second peak was related to developing a DCM at the metalimnion during the stratification reaching maximum concentrations in September-October. average. The above and below the average Chl-a concentration regularly alternated (three above, two below, three above, three below, two above, four below, two above, two below). Annually, Chl-a concentration displayed a regular dynamic with two peaks ( Figure  3). The first peak was associated with the winter (January to March) mixing period when diatoms bloom throughout the water column. The second peak was related to developing a DCM at the metalimnion during the stratification reaching maximum concentrations in September-October.

Euphotic Zone (ZEU)
There was interannual variability in the euphotic zone (ZEU) of Lake Alchichica from 1998 to 2018 ( Figure 4). A recurrent annual pattern evidences the turbid and the clear water phases. There was a linearly decreasing trend in the ZEU from 1998 to 2018. The ZEU annual average was 31 ± 6 m within 27 ± 6 m and 37 ± 8 m.

Euphotic Zone (Z EU )
There was interannual variability in the euphotic zone (Z EU ) of Lake Alchichica from 1998 to 2018 ( Figure 4). A recurrent annual pattern evidences the turbid and the clear water phases. There was a linearly decreasing trend in the Z EU from 1998 to 2018. The Z EU annual average was 31 ± 6 m within 27 ± 6 m and 37 ± 8 m.  The ZEU fluctuated interannually ( Figure 5). For 29%, the ZEU was above the annual average, while for 71%, it was below the annual average. The ZEU displayed a regular dynamic with two phases annually ( Figure 5). The first is the turbid water phase with ZEU below the annual average (31 m) during the mixing and early stratification (January to June). For the second one, the clear water phase with ZEU was above the annual average during the well-established and late stratification (July to December). The Z EU fluctuated interannually ( Figure 5). For 29%, the Z EU was above the annual average, while for 71%, it was below the annual average. The Z EU displayed a regular dynamic with two phases annually ( Figure 5). The first is the turbid water phase with Z EU below the annual average (31 m) during the mixing and early stratification (January to June). For the second one, the clear water phase with Z EU was above the annual average during the well-established and late stratification (July to December).

Mixed Layer Depth (Z MIX )
The mixed layer depth (Z MIX ) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 6). A recurrent annual pattern evidences the water column mixing and stratification periods. The Z MIX annual average was 29 ± 10 m within 24 ± 16 m and 34 ± 21 m. There was a linearly decreasing trend in the Z MIX from 1998 to 2018.
The ZEU fluctuated interannually ( Figure 5). For 29%, the ZEU was above the annual average, while for 71%, it was below the annual average. The ZEU displayed a regular dynamic with two phases annually ( Figure 5). The first is the turbid water phase with ZEU below the annual average (31 m) during the mixing and early stratification (January to June). For the second one, the clear water phase with ZEU was above the annual average during the well-established and late stratification (July to December).

Mixed Layer Depth (ZMIX)
The mixed layer depth (ZMIX) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 6). A recurrent annual pattern evidences the water column mixing and stratification periods. The ZMIX annual average was 29 ± 10 m within 24 ± 16 m and 34 ± 21 m. There was a linearly decreasing trend in the ZMIX from 1998 to 2018. The ZMIX fluctuated interannually ( Figure 7). For 52% of the period, the ZMIX was above the annual average, while 48% was below the annual average. Annually, the ZMIX displayed a regular dynamic with two seasons. The first is when ZMIX is maximum (i.e., the whole water column) during the winter mixing period (January to March). The second one corresponded to the extended stratification season (April to December). The ZMIX decreased rapidly in April at the onset of the stratification and from then on increased gradually along with the well-established and late stratification. The ZMIX was above the annual average from January to March and November to December and below the average from April to October (Figure 7). The Z MIX fluctuated interannually ( Figure 7). For 52% of the period, the Z MIX was above the annual average, while 48% was below the annual average. Annually, the Z MIX displayed a regular dynamic with two seasons. The first is when Z MIX is maximum (i.e., the whole water column) during the winter mixing period (January to March). The second one corresponded to the extended stratification season (April to December). The Z MIX decreased rapidly in April at the onset of the stratification and from then on increased gradually along with the well-established and late stratification. The Z MIX was above the annual average from January to March and November to December and below the average from April to October (Figure 7).

Temperature of the Mixed Layer (T MIX )
The mixed layer temperature (T MIX ) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 8). A recurrent annual pattern evidences the cold/dry and the warm/rainy seasons associated with the tropical climate. The T MIX annual average was 17.4 ± 1.8 • C within a range of 16.9 ± 1.7 • C and 18.2 ± 1.8 • C. There was an increasing linear trend in the T MIX from 1998 to 2018. displayed a regular dynamic with two seasons. The first is when ZMIX is maximum (i.e., the whole water column) during the winter mixing period (January to March). The second one corresponded to the extended stratification season (April to December). The ZMIX decreased rapidly in April at the onset of the stratification and from then on increased gradually along with the well-established and late stratification. The ZMIX was above the annual average from January to March and November to December and below the average from April to October (Figure 7).

Temperature of the Mixed Layer (TMIX)
The mixed layer temperature (TMIX) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 8). A recurrent annual pattern evidences the cold/dry and the warm/rainy seasons associated with the tropical climate. The TMIX annual average was 17.4 ± 1.8 °C within a range of 16.9 ± 1.7 °C and 18.2 ± 1.8 °C. There was an increasing linear trend in the TMIX from 1998 to 2018. The TMIX fluctuated interannually ( Figure 9). For 48% of the period, the TMIX was above the annual average, while 52% was below the annual average. The TMIX displayed a regular dynamic with two periods (Figure 9). The first was when TMIX was colder (below the average) during the winter and late stratification (January-March, and November-December). The second one corresponded to the well-established and late stratification season (April to October) when TMIX was warmer (above the average). The T MIX fluctuated interannually ( Figure 9). For 48% of the period, the T MIX was above the annual average, while 52% was below the annual average. The T MIX displayed a regular dynamic with two periods (Figure 9). The first was when T MIX was colder (below the average) during the winter and late stratification (January-March, and November-December). The second one corresponded to the well-established and late stratification season (April to October) when T MIX was warmer (above the average).

Dissolved Oxygen of the Mixed Layer (DO MIX )
The dissolved oxygen concentration of the mixed layer (DO MIX ) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 10). A recurrent annual pattern evidences the water column mixing and stratification periods. The DO MIX annual average was 6.3 ± 1.0 mg L −1 within 5.6 ± 1.1 mg L −1 and 6.9 ± 0.9 mg L −1 . There was an increasing linear trend in the DO MIX from 1998 to 2018.
The TMIX fluctuated interannually (Figure 9). For 48% of the period, the TMIX was above the annual average, while 52% was below the annual average. The TMIX displayed a regular dynamic with two periods (Figure 9). The first was when TMIX was colder (below the average) during the winter and late stratification (January-March, and November-December). The second one corresponded to the well-established and late stratification season (April to October) when TMIX was warmer (above the average).

Dissolved Oxygen of the Mixed Layer (DOMIX)
The dissolved oxygen concentration of the mixed layer (DOMIX) displayed interannual variability in Lake Alchichica from 1998 to 2018 ( Figure 10). A recurrent annual pattern evidences the water column mixing and stratification periods. The DOMIX annual average was 6.3 ± 1.0 mg L −1 within 5.6 ± 1.1 mg L −1 and 6.9 ± 0.9 mg L −1 . There was an increasing linear trend in the DOMIX from 1998 to 2018. The DOMIX fluctuated interannually ( Figure 11). For 57% of the period, the DOMIX was above the annual average, while 43% was below the annual average. The DOMIX displayed a regular seasonal dynamic with two periods (Figure 11). The first was when DOMIX was lower (below the annual average) during the winter and late stratification (January-March and November-December). The second one corresponded to the well-established and late stratification season (April to October) when DOMIX was higher (above the annual average). The reduction in DOMIX takes place from the late stratification (through the erosion of the thermocline) to the mixing period when the anoxic hypolimnetic water mixes with the aerobic epilimnetic waters. The DO MIX fluctuated interannually ( Figure 11). For 57% of the period, the DO MIX was above the annual average, while 43% was below the annual average. The DO MIX displayed a regular seasonal dynamic with two periods (Figure 11). The first was when DO MIX was lower (below the annual average) during the winter and late stratification (January-March and November-December). The second one corresponded to the well-established and late stratification season (April to October) when DO MIX was higher (above the annual average). The reduction in DO MIX takes place from the late stratification (through the erosion of the thermocline) to the mixing period when the anoxic hypolimnetic water mixes with the aerobic epilimnetic waters.

Dissolved Inorganic Nitrogen of the Mixed Layer (DIN MIX )
The dissolved inorganic nitrogen concentration in the mixed layer (DIN MIX ) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure 12). There was no clear annual pattern of DIN MIX associated with the water column mixing and stratification periods. The DIN MIX annual average was 2.7 ± 2.0 µM within a range of 1.1 ± 0.5 µM and 5.1 ± 3.4 µM. There was an increasing linear trend in the DIN MIX from 1999 to 2018. lower (below the annual average) during the winter and late stratification (January-March and November-December). The second one corresponded to the well-established and late stratification season (April to October) when DOMIX was higher (above the annual average). The reduction in DOMIX takes place from the late stratification (through the erosion of the thermocline) to the mixing period when the anoxic hypolimnetic water mixes with the aerobic epilimnetic waters.

Dissolved Inorganic Nitrogen of the Mixed Layer (DINMIX)
The dissolved inorganic nitrogen concentration in the mixed layer (DINMIX) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure  12). There was no clear annual pattern of DINMIX associated with the water column mixing and stratification periods. The DINMIX annual average was 2.7 ± 2.0 µM within a range of 1.1 ± 0.5 µM and 5.1 ± 3.4 µM. There was an increasing linear trend in the DINMIX from 1999 to 2018. The DINMIX fluctuated interannually ( Figure 13). For 45% of the period, the DINMIX was above the annual average, while 55% was below the annual average. The DINMIX displayed a regular dynamic with two periods (Figure 13). The first was when DIN was lower (below the annual average) at the onset and along most of the stratification (April, June to October, and December). The second one corresponded primarily to the mixing season (January to March), May, and November, when DIN was higher (above the annual average). The DIN MIX fluctuated interannually ( Figure 13). For 45% of the period, the DIN MIX was above the annual average, while 55% was below the annual average. The DIN MIX displayed a regular dynamic with two periods (Figure 13). The first was when DIN was lower (below the annual average) at the onset and along most of the stratification (April, June to October, and December). The second one corresponded primarily to the mixing season (January to March), May, and November, when DIN was higher (above the annual average).

Soluble Reactive Phosphorus of the Mixed Layer (SRP MIX )
The soluble reactive phosphorous concentration in the mixed layer (SRP MIX ) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure 14). There was no clear annual pattern of the SRP MIX associated with the water column mixing and stratification periods. The SRP MIX annual average was 0.27 ± 0.16 µM within a range of 0.15 ± 0.11 µM and 0.53 ± 0.24 µM.
The DINMIX fluctuated interannually ( Figure 13). For 45% of the period, the DINMIX was above the annual average, while 55% was below the annual average. The DINMIX displayed a regular dynamic with two periods (Figure 13). The first was when DIN was lower (below the annual average) at the onset and along most of the stratification (April, June to October, and December). The second one corresponded primarily to the mixing season (January to March), May, and November, when DIN was higher (above the annual average).

Soluble Reactive Phosphorus of the Mixed Layer (SRPMIX)
The soluble reactive phosphorous concentration in the mixed layer (SRPMIX) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure  14). There was no clear annual pattern of the SRPMIX associated with the water column mixing and stratification periods. The SRPMIX annual average was 0.27 ± 0.16 µM within a range of 0.15 ± 0.11 µM and 0.53 ± 0.24 µM. The SRPMIX fluctuated widely interannually ( Figure 15). For 40% of the period, the SRPMIX was above the annual average, and for 60% was below the annual average. Seasonally, the SRPMIX displayed a regular dynamic with two periods (Figure 15). The first was when SRPMIX was lower (below the annual average), and most of the stratification (June to December). The second one corresponded to the mixing season and early stratification (January to May) when SRPMIX was higher (above the annual average).  The SRP MIX fluctuated widely interannually ( Figure 15). For 40% of the period, the SRP MIX was above the annual average, and for 60% was below the annual average. Seasonally, the SRP MIX displayed a regular dynamic with two periods (Figure 15). The first was when SRP MIX was lower (below the annual average), and most of the stratification (June to December). The second one corresponded to the mixing season and early stratification (January to May) when SRP MIX was higher (above the annual average).

Soluble Reactive Silica of the Mixed Layer (SRSi MIX )
The soluble reactive silica concentration in the mixed layer (SRSi MIX ) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure 16). There was not a clear annual pattern of the SRSi MIX associated with the water column mixing and stratification periods. The SRSi MIX annual average was 4.6 ± 3.6 µM within a range of 1.9 ± 1.1 µM and 11.6 ± 5.2 µM.
The SRPMIX fluctuated widely interannually ( Figure 15). For 40% of the period, the SRPMIX was above the annual average, and for 60% was below the annual average. Seasonally, the SRPMIX displayed a regular dynamic with two periods (Figure 15). The first was when SRPMIX was lower (below the annual average), and most of the stratification (June to December). The second one corresponded to the mixing season and early stratification (January to May) when SRPMIX was higher (above the annual average).

Soluble Reactive Silica of the Mixed Layer (SRSiMIX)
The soluble reactive silica concentration in the mixed layer (SRSiMIX) displayed substantial interannual variability in Lake Alchichica from 1999 to 2018 ( Figure 16). There was not a clear annual pattern of the SRSiMIX associated with the water column mixing and stratification periods. The SRSiMIX annual average was 4.6 ± 3.6 µM within a range of 1.9 ± 1.1 µM and 11.6 ± 5.2 µM. The SRSiMIX fluctuated extensively interannually ( Figure 17). For 30% of the period, the SRSiMIX was above the annual average, and for 70% was below the annual average. The SRSiMIX displayed irregular dynamic alternating values below and above the annual average ( Figure 17). One-third of the year, SRSiMIX was below the annual average (February, May-June, and September). For the rest of the year (January, March-April, July-August, and October-December), the SRSiMIX was above the annual average.  The SRSi MIX fluctuated extensively interannually ( Figure 17). For 30% of the period, the SRSi MIX was above the annual average, and for 70% was below the annual average. The SRSi MIX displayed irregular dynamic alternating values below and above the annual average ( Figure 17). One-third of the year, SRSi MIX was below the annual average (February, May-June, and September). For the rest of the year (January, March-April, July-August, and October-December), the SRSi MIX was above the annual average.

Wavelet Analyses
Wavelet analysis was used to decompose the Chl-a variability from 1998 to 2018 as a function of period and time. We examined if the continuous wavelet transforms could successfully identify variability in the annual scale. The local wavelet power spectrum measures the variance distribution of the time series according to time and periodicity; high variability is represented by red color, whereas blue indicates weak variability. The wavelet spectrum showed that Chl-a ( Figure 18) displayed a consistent 1-year periodicity from 1998 until 2012. The strength of the annual periodicity changed over time, suggesting a non-stationary behavior of the time series. There was a weakening of the annual periodicity from 2008 to 2012 until it became statistically non-significant from 2013 to 2015. Moreover, a consistent and statistically significant periodicity between~4 and 5 years in Chl-a was detected. A semi-annual periodicity was significant for a short time, around 2016, but did not emerge as significant in the global wavelet spectrum. Wavelet analysis of the Multivariate "El Niño Southern Oscillation" (ENSO) Index (MEI) time series (Figure 18) showed a 2-year and a four to five-year periodicity. The SRSiMIX fluctuated extensively interannually (Figure 17). For 30% of the period, the SRSiMIX was above the annual average, and for 70% was below the annual average. The SRSiMIX displayed irregular dynamic alternating values below and above the annual average ( Figure 17). One-third of the year, SRSiMIX was below the annual average (February, May-June, and September). For the rest of the year (January, March-April, July-August, and October-December), the SRSiMIX was above the annual average.

Wavelet Analyses
Wavelet analysis was used to decompose the Chl-a variability from 1998 to 2018 as a function of period and time. We examined if the continuous wavelet transforms could successfully identify variability in the annual scale. The local wavelet power spectrum measures the variance distribution of the time series according to time and periodicity; high variability is represented by red color, whereas blue indicates weak variability. The wavelet spectrum showed that Chl-a ( Figure 18  Chl-a was detected. A semi-annual periodicity was significant for a short time, around 2016, but did not emerge as significant in the global wavelet spectrum. Wavelet analysis of the Multivariate "El Niño Southern Oscillation" (ENSO) Index (MEI) time series ( Figure  18) showed a 2-year and a four to five-year periodicity. However, the former 2-year periodicity showed shifts between 2002 and 2004 and between 2012 and 2015, whereas the latter's (~4-5 years) period remained stable for the duration of the time series. This shift in periodicity coincides with the shifts in periodicity displayed by the Chl-a wavelet spectrum (Figure 9).Wavelet analysis of the ZEU, ZMIX, TMIX, and DOMIX showed stable 1-year periodicities for the whole duration of the time series ( Figure 19). This steady behavior means that the system's complexity has not deteriorated so that the seasonality of this process has not perceptibly changed over time. In addition, a 4-year periodicity was significant for the ZMIX time series. Figure 18. Continuous wavelet power spectra showing the periodicity (left) and time-averaged (above global) wavelet power spectra (right) of (a) Chl-a concentration integrated per unit area in Lake Alchichica. Continuous wavelet power spectra showing the periodicity (left) and time-averaged (above global) wavelet power spectra (right) of (b) the Multivariate ENSO Index (MEI) in Lake Alchichica.
Finally, the wavelet power spectrum for nutrients (DINMIX = dissolved inorganic nitrogen, SRPMIX = soluble reactive phosphorus, SRSiMIX = soluble reactive silica) also showed statistically significant periodic structures by the annual (1-year) cycle and, in addition, a 2-year and a ~4-5 years statistically significant periodicities were also characterized for the time series ( Figure 20). All significant scales in the wavelet spectrum were clustered around the 1-year period. However, the power of the annual scale was much lower from 2002 to 2004, except for the SRPMIX time series. Four-year periodicities were statistically significant for most of the nutrient time series. However, the DINMIX time series 2-year periodicity became  Figure 19). This steady behavior means that the system's complexity has not deteriorated so that the seasonality of this process has not perceptibly changed over time. In addition, a 4-year periodicity was significant for the Z MIX time series.  Finally, the wavelet power spectrum for nutrients (DIN MIX = dissolved inorganic nitrogen, SRP MIX = soluble reactive phosphorus, SRSi MIX = soluble reactive silica) also showed statistically significant periodic structures by the annual (1-year) cycle and, in addition, a 2-year and a~4-5 years statistically significant periodicities were also characterized for the time series (Figure 20).

Discussion
Winder and Cloern [20] found the most common pattern reported in temperate and subtropical aquatic ecosystems, one phytoplankton peak per year or the "canonical spring-bloom pattern" [5,36], to be less frequent than expected (48%), and although the climate-driven annual cycle is an important driver, other factors can mask it leading to other periodicities or even the lack of periodicity [14]. Melack [1] found that most tropical lakes exhibit seasonal fluctuations associated with rainfall, river discharge, or vertical mixing, corresponding to the hydrological (water input-output) or hydrographic (water column and circulation) features of Talling [10]. Lake Alchichica has no river discharge and

Discussion
Winder and Cloern [20] found the most common pattern reported in temperate and subtropical aquatic ecosystems, one phytoplankton peak per year or the "canonical spring-bloom pattern" [5,36], to be less frequent than expected (48%), and although the climatedriven annual cycle is an important driver, other factors can mask it leading to other periodicities or even the lack of periodicity [14]. Melack [1] found that most tropical lakes exhibit seasonal fluctuations associated with rainfall, river discharge, or vertical mixing, corresponding to the hydrological (water input-output) or hydrographic (water column and circulation) features of Talling [10]. Lake Alchichica has no river discharge and lies in a semi-desertic area; then, the seasonality is driven by hydrographic features controlled by its thermal stratification and mixing pattern.
Deep tropical lakes are warm-monomictic, showing great regularity in seasonal mixing, which typically coincides with the hemispheric winter (e.g., [37]). Previous studies on Lake Alchichica suggested the latter (e.g., [22,38,39]) and herein confirmed during the 21 annual cycles studied. Lake Alchichica showed a remarkable regularity on the hydrodynamical pattern (warm monomixis), displaying a predictable annual, complete mixing period in winter (January to March) and remained stratified for the remainder of the year. Adame et al. [3] subdivided the stratification period into three stages: (a) early stratification (April to June) when the thermocline was close to the surface and poorly defined, (b) well-established stratification (July to September) when the water column displays a clear layering (epi, meta, and hypolimnion), and (c) late stratification (October to December) when the thermocline deepened and the metalimnion thinned.
The recurrent annual thermal stratification and mixing pattern of Lake Alchichica and its main limnological features could be summarised as follow (Table 1). Z MIX is at maximum (from surface down to the bottom) during the mixing season, and the water column is cold and well-oxygenated. Nutrients (DIN, SRP, and SRSi) reach high concentrations along the water column, promoting the winter diatom bloom development, resulting in high Chl-a concentration and the lowest values of Z EU . The mixing season is the turbid water phase associated with high biogenic turbidity. The winter bloom comprises the centric diatom Cyclotella alchichicana [3,26].
The stratification period corresponds to the clear water phase with high values of Z EU related to the presence of low Chl-a concentration in the Z MIX . The Z MIX remains warm and well-oxygenated. Through the extended stratification season, Z MIX and Z EU broaden.
While stratified, nutrients are rapidly depleted in the Z MIX and remain so for the remaining of the stratification, explaining the low Chla concentration in the Z MIX and the high values of Z EU .
The large size of C. alchichicana makes it inedible for zooplankton [40,41]. Then most of the diatom bloom is not consumed and incorporated in the food web but instead exported below the thermocline, down to the bottom [42]. A large amount of organic matter deposited into the sediments is partially remineralized by consuming the available dissolved oxygen leading to hypolimnetic anoxia [3], and partially buried and sequestered into the sediments [43]. During the stratification period, the organic matter remineralization explains the nutrient accumulation in the hypolimnion. Along the stratification period, nutrients are depleted in the epilimnion and accumulated in the hypolimnion.
Approximately from September to November, a DCM develops at the metalimnion. The Chl-a concentration at the DCM reaches high values, similar to those recorded in the winter diatom bloom. The DCM is mainly composed of C. alchichicana, which is exported below the thermocline down to the sediments [42]. In this way, two Chl-a peaks are recognized in the annual cycle: a diatom winter blooms along the water column during the mixing season. The second is a DCM at the metalimnion in the stratification.
As previously mentioned, Lake Alchichica showed not one but two phytoplankton peaks per year ( Figure 21): the winter bloom (Chl-a: 23 ± 5 mg m −2 ) and the DCM (Chl-a: 28 ± 12 mg m −2 ). The average Chl-a concentration of the mixing period (January to March) was higher than in the DCM in five years (24%), while the Chl-a concentration in the DCM was higher than in the mixing season in sixteen years (76%). The winter diatom bloom duration fluctuated between 1 and 4 months (January-April), with Chl-a concentrations above the 21-year average. January (11, 52. exported below the thermocline, down to the bottom [42]. A large amount of organic matter deposited into the sediments is partially remineralized by consuming the available dissolved oxygen leading to hypolimnetic anoxia [3], and partially buried and sequestered into the sediments [43]. During the stratification period, the organic matter remineralization explains the nutrient accumulation in the hypolimnion. Along the stratification period, nutrients are depleted in the epilimnion and accumulated in the hypolimnion. Approximately from September to November, a DCM develops at the metalimnion. The Chl-a concentration at the DCM reaches high values, similar to those recorded in the winter diatom bloom. The DCM is mainly composed of C. alchichicana, which is exported below the thermocline down to the sediments [42]. In this way, two Chl-a peaks are recognized in the annual cycle: a diatom winter blooms along the water column during the mixing season. The second is a DCM at the metalimnion in the stratification. As previously mentioned, Lake Alchichica showed not one but two phytoplankton peaks per year ( Figure 21): the winter bloom (Chl-a: 23 ± 5 mg m −2 ) and the DCM (Chl-a: 28 ± 12 mg m −2 ). The average Chl-a concentration of the mixing period (January to March) was higher than in the DCM in five years (24%), while the Chl-a concentration in the DCM was higher than in the mixing season in sixteen years (76%). The winter diatom bloom duration fluctuated between 1 and 4 months (January-April), with Chl-a concentrations above the 21-year average. January (11, 52.  Lewis [44] compared the phytoplankton biomass in temperate and tropical lakes by analyzing the ratio of annual and seasonal mean to maximum biomass. Tropical lakes displayed a higher ratio (0.45) of annual mean to maximum than temperate lakes (0.36). However, the ratio of seasonal mean to maximum biomass (~0.7) was similar in both tropical and temperate lakes. The latter means the intensity of the annual bloom in both tropical and temperate lakes are similar, but tropical lakes display higher minimum Chl-a con- Lewis [44] compared the phytoplankton biomass in temperate and tropical lakes by analyzing the ratio of annual and seasonal mean to maximum biomass. Tropical lakes displayed a higher ratio (0.45) of annual mean to maximum than temperate lakes (0.36). However, the ratio of seasonal mean to maximum biomass (~0.7) was similar in both tropical and temperate lakes. The latter means the intensity of the annual bloom in both tropical and temperate lakes are similar, but tropical lakes display higher minimum Chl-a concentrations throughout the year. The annual ratio of mean to maximum Chl-a in Lake Alchichica ranged between 0.36 and 0.85 (0.68 ± 0.13), which is higher than those reported by Lewis [44] for tropical lakes (0.39 to 0.58) but similar to the seasonal ratio of mean to maximum Chl-a with 0.65 to 0.80. The discrepancy could be that Lewis [44] calculated the seasonal ratio by considering the 3-month interval with the highest biomass, most likely associated with the winter (mixing) bloom. The seasonal ratio of mean to maximum Chl-a calculated for the winter diatom bloom (3 months, January to March) in Lake Alchichica is 0.62, slightly lower than the seasonal ratio of Lewis [44] with 0.65 to 0.80.
The seasonal ratio of mean to maximum Chl-a calculated for the DCM (4 months, August to November) is 0.36, much lower than the seasonal ratio of Lewis [44]  The annual variation in phytoplankton biomass in temperate lakes is greater than in tropical lakes [44]. Kalff and Watson [45] found that Lake Naivasha had low and seasonal variable (10×) phytoplankton biomass. Differently, the coefficient of variation in Lake Alchichica annual phytoplankton biomass averages 27 ± 14%, and low (1.3×) seasonal variation in biomass (i.e., mixing: 23 ± 5 mg m −2 , early stratification: 19 ± 4 mg m −2 , well-established stratification: 25 ± 10 mg m −2 , and late stratification: 25 ± 10 mg m −2 ). There were also interannual differences in the max:min biomass ratio with an average of 2.4 ± 1.0 (1.5 to 6.1). The inter-annual variability of Chl-a was high -.7×-(annual mean: 18 to 30 mg m −2 ), similarly to Lake Kivu with 1.9× [12].
Our analyses showed that continuous wavelet transforms successfully detect the phenology of the Chl-a concentration and associated limnological parameters over time. Throughout the time series, periodicities associated with the annual cycle represented the dominant scale of variation. As mentioned, the recurrent and predictable annual cyclicity is explained mainly by the warm-monomictic thermal pattern of Lake Alchichica, with the water column mixing being the most important driving force explaining the seasonality.
However, the Chl-a continuous wavelet transforms identified a temporary decrease in the strength of the annual periodicity. While the annual periodicity was the most important scale of variation for the Chl-a during 1998-2010, it was significantly weakened from 2011 to 2015, when it became no longer statistically significant. In this period, there were slight decreases (Z EU , Z MIX ) and increases (T MIX , DO MIX ) in the limnological variables associated with decreases in Chl-a concentration. Interestingly, this shift in periodicity in the Chl-a power spectrum between 2011 and 2015 coincided with intense La Niña events (2011-2012) and with the "Godzilla" El Niño event (2015).
Tropical lakes are particularly sensitive to climate change [46,47]. A previous study [40]) reported that the 1998-1999 El Niño resulted in higher epilimnetic temperatures and lower Chl-a concentrations in Lake Alchichica. Recently, Alcocer et al. [40] found that the air temperature in the Lake Alchichica region increased from 1966 to 2018. This variability was correlated positively with the Oceanic Niño Index (ONI), the Pacific Decadal Oscillation index (PDO), and the Atlantic Multidecadal Oscillation index (AMO), and negatively with the Southern Oscillation Index (SOI). T MIX temperature displayed an increasing trend associated with the air temperature increase, while Chl-a concentration displayed a decreasing trend (Table 2). Higher temperatures associated with climate change resulted in stronger thermal stratifications, lower Z MIX, and lower nutrient concentration availability. DIN MIX and SRP MIX displayed decreasing trends; contradictorily, SRSi MIX showed an increasing trend related to the decreasing trend in Chl-a, indicating lower biomass of diatoms and consequently less SRSi uptake. Regarding the 2-year and 4-year periodicities revealed through the continuous wavelet transforms analyses, the biennial cycle in the Chl-a concentration was already reported in Adame et al. [3]. Adame et al. [3] reported higher Chl-a concentrations during the winter diatom blooms of even years while lower in odd years. Higher silicate concentration and Si/N ratio in the hypolimnion during the previous late stratification period led to higher diatom blooms in the following mixing period. Herein, we confirmed the biennial (2-year) cyclicity and provided a plausible explanation of the role silicate and the Si/N ratio play in controlling the phytoplankton biomass in Lake Alchichica.
As mentioned before, the large size of C. alchichicana makes it inedible for zooplankton. Then most of the biomass is exported below the thermocline down to the sediments [27]. C. alchichicana contributes 98% to the silica flux down to the sediments in Lake Alchichica with 147 g m −2 year −1 [42]. Oseguera et al. [48] found the particles captured in sediment traps deployed in Lake Alchichica evaluating exportation into the deep sediments were biogenic, primarily composed of cells of C. alchichicana.
Alcocer et al. [43] found high preservation of biogenic carbon-and associated diatom silica-in the deep sediments of Lake Alchichica, suggesting a net carbon and silica lost out of the system. Lake Alchichica, with a low silica concentration availability and phytoplankton biomass dominated by large diatoms, deep, thorough mixing events are crucial for silica mineralization and recycling instead of sequestering. Mineralization and mixing events also replenish the lake with nitrogen which often limits phytoplankton in tropical lakes [49], as in Lake Alchichica [25,50], confirmed by the recurrent annual bloom of N. aff. spumigena [28], a nitrogen-fixing cyanobacteria [51].
Climate variability leads to warmer temperatures, which increase the length and strength of thermal stratification (e.g., [52]), preventing nutrient mineralization and recycling. Opposite, colder temperatures favor complete water column mixing and reoxygenation, favoring nutrient mineralization and recycling, which brought an increased supply of phosphorus and, more importantly, nitrogen which later controlled phytoplankton growth in the Z EU .
Moreover, there was a~4-5 years periodicity evidenced in nutrient concentration. We think these 2-year and~4-5 years periodicities could be associated with ENSO events encompassing stronger winds and then more profound, more complete mixing seasons providing higher nutrient concentrations stimulating larger blooms as suggested by Alcocer et al. [40].

Conclusions
This work confirmed that wavelet analysis could help us interpret multi-scale, nonstationary time-series data to reveal features that could not otherwise be seen. By using a long-term (1998-2018) limnological database of the tropical, saline Lake Alchichica, we (a) verified the regularity of the annual cycle of the phytoplankton Chl-a associated with the warm-monomictic thermal-mixing-pattern, (b) revealed the phenology of the annual Chl-a cycle showing large inter-annual variability, (c) identified the presence of a biennial cycle most likely controlled through SRSi availability and Si/N ratio associated to climate variability influencing the length and strength of the annual thermal stratification, and (d) recognize other cycles with different periodicities associated with external forcing agents (e.g., ENSO-El Niño Southern Oscillation). The water quality variables were regular along the 21 years, displaying a recurrent annual pattern. Nutrient concentration, as trophic indicators, revealed the same cyclicity (1-year, 2-year, and~4-5 years) as Chl-a, suggesting the external forcing agents lead to deeper, more complete mixing events increasing nutrient availability, which in turn stimulate larger Chl-a concentrations.