Fúquene Lake is a large (30 km²) and shallow (average depth 2 m, maximum depth 6 m), Andean lake (2,560 m a.s.l.) located on the eastern Cordillera of the Colombian Andes (5°27′55 ′′N, 75°46′19′′ W) (Figure 1). The Ubaté, Lenguazaque and Suta rivers are the main inflows and the outflow to the north forms the Suárez River, a key tributary of the Magdalena River. Fúquene Lake is an important provider of natural habitats for native aquatic birds, fish, and macrophytes and offers many benefits to local inhabitants including water and food (Vidal and Perez, 2007).

Fúquene Lake was once regarded as one of the largest natural Andean shallow lakes in Colombia (approx. area of 112 km²; Franco et al., 2011). Archaeological records and pollen data from the lake catchment area have shown that humans occupied and modified the landscape for agriculture well before 1500 CE (Van Geel and Van der Hammen, 1973; Langebaek, 1995). After the Spanish conquest at 1500 CE there were major increases in population densities at nearby settlements (Langebaek, 1995) marked by large-scale land-use shifts from small-cultivated areas into large grazing areas (Etter et al., 2008). Europeans introduced several non-native grazing animals such as cattle, goats, sheep, and horses, which quickly transformed the landscape from natural montane forests/cultivated areas to large extensions of pasture lands (Etter et al., 2008). Since the late 1800s Fúquene Lake was further subjected to several draining episodes which aimed to increase land for agriculture and livestock. By the mid-1930s the lake had lost almost 70% of its original surface area (Figure 1; Table 1; Vidal and Perez, 2007). More recently, in 1987, the lake hydrology was further affected through the construction of a canal along the lake shorelines which was built to stop further expansion of agricultural lands (Franco et al., 2011). This engineering work impeded lateral hydrological connectivity of the lake, diverted most of the inflow waters, and increased the water retention time of the lake, which are thought to have stimulated lake productivity and the spread of invasive species (Table 1; Vidal et al., 2015). The lake is currently in a hypertrophic state (total phosphorus [TP] > 0.40 mg/L) with well-established populations of Eichhornia and Egeria (Franco et al., 2011; Figure 1). The date of introduction of Eichhornia, which is believed to be a native plant to the lowland Amazon and lower Orinoco basins (Barrett and Forno, 1982), is not fully clear but locals affirm that the plant was transplanted into the lake sometime around the mid-1950s. In contrast, E. densa, native to the southern Hemisphere (Yarrow et al., 2009), was deliberately introduced in 1990 as a phytoremediation strategy (Duque and Donato, unpublished) to improve the poorly oxygenated (dissolved oxygen of 4.41 ppm) eutrophic (total phosphorus-TP of 1.73 mg/L) waters at the time (Table 1; Donato et al., 1987).

We selected four sampling transects of approx. 1 km long across the lake to characterize the present-day macrophyte species composition and environmental conditions (Figure 1). Each transect was defined according to the dominant macrophyte type and designated with a code (in parentheses): (i) E. crassipes (WH1); (ii) E. densa (ED1); (iii) mix assemblage of macrophyte species (M1); and (iv) open waters (OW1) (Figure 1). At each transect we collected macrophyte data at approximately 50 m intervals (n = 18 sampling points), during five fieldwork sessions (December 2015–April 2017). Surveys were undertaken from a boat and at each sampling point we used a combination of bathyscope and double-headed rake sampler, and macrophyte species cover percentages within a 1 m² area were recorded. Due to difficulties in accessing inner areas of the dense floating mats of Eichhornia, macrophyte-sampling points at transect WH1 were assessed from the margins. Concurrent with the macrophyte surveys, we collected three measurements of water depth, secchi depth, water temperature, conductivity, salinity, total dissolved solids (TDS) and dissolved oxygen (DO) at the beginning, middle, and end of each transect. For transect WH1, we removed or displaced a few plant individuals of Eichhornia at each environmental sampling point and waited for a few minutes to let hydrological conditions settle back before measuring the limnological variables. Physical and chemical parameters were measured at the surface water and at 0.5 m depth using an YSI® DO and conductivity meter.

We used a corer of 10 cm diameter (Aquatic Instruments Inc.) to collect a large volume of sediment whilst preserving undisturbed surface sediments (Patmore et al., 2014). We retrieved two short sediment cores with lengths of 100 cm (LFUQ-M1) and 50 cm (LFUQ-WH1) from two of the contemporary sampling transects (M1, and WH1, respectively; Figure 1C). The core LFUQ-M1 was retrieved at the end of transect M1 (5°27′13.86′′N, 73°44′59.33′′W) at a water depth of 100 cm and core LFUQ-WH1 was collected at the mid-point of the transect at a water depth of 90 cm (5°27′45.35′′N, 73°45′48.61′′W). Each core was extruded in the field at 1 cm intervals and changes in lithostratigraphy were recorded prior to extrusion.

We dated the longer core LFUQ-M1 by combining ²¹⁰Pb and ¹⁴C. Radionuclide measurements of ²¹⁰Pb (half-life 22.3 years), ¹³⁷Cs and ²⁴¹Am were performed on the most recent samples (0–20 cm; Appleby et al., 1986). Dates were ascribed using the Constant Rate of Supply (CRS) model (Appleby and Oldfield, 1978). ¹⁴C analysis on plant macrofossil remains via accelerator mass spectrometry analysis (AMS; Beta Analytic Inc.) was used to date older sediments at depths of 30, 75, and 98 cm. A combined age-depth model was calculated using the CLAM 2.2 software in R (Blaauw, 2010; R Development Core Team, 2016) to ascribe dates below the ²¹⁰Pb age model. Radiocarbon ages were calibrated using the Northern Hemisphere calibration curve and an age model used linear loess interpolation (lower AIC and no time reversals). The final age model that best matched the environmental history of the lake resulted from using the maximum error standard ²¹⁰Pb dates for the top 15 cm intervals and the minimum error standard dates derived from the CLAM model for 15–100 cm sediment intervals (Figure 3; Table 2).

Geochemical element concentrations were measured on 1 cm-thick discrete samples using X-Ray Fluorescence (XRF). For the cores LFUQ-M1 sampling resolution was at 2-cm intervals (0–20 cm) and at 4-cm intervals from 20 to 100 cm. For the core LFUQB1 we analyzed contiguous sediment intervals (0–20 cm) and every 3-cm between 20 and 50 cm. The XRF analyses were measured in the laboratory using a portable XRF analyzer spectrometer, XMET 7500. Dry sediment samples were ground and homogenized in the laboratory using a mortar and pestle. Approximately 3 g of sediment sample were used and covered with Chemplex thin-film sample support. Mean values of each element were determined and plotted using two readings per sample. Calcium (Ca), Potassium (K), Iron (Fe), Manganese (Mn), Rubidium (Rb), and Titanium (Ti) results were selected for this study. Given that Ti is an unambiguous indicator of allochthonous inputs from the catchment, we used this element for normalization (Cohen, 2003). We calculated four different but complementary index ratios according to Davies et al. (2015) to investigate changes in: weathering regimes (Rb/K), physical erosion of the lake catchment area (K/Ti), lake level fluctuation (Ca/Ti) and lake reduction conditions ([Fe/Mn]/Ti). We also used loss-on-ignition (LOI) (Dean, 1974) to determine the organic fraction of contiguous sediment intervals (0–50 cm) in the cores LFUQ-M1 and LFUQ-WH1, and of 2-cm intervals from 50 to 100 cm in core LFUQ-M1.

As funds were not available for dating core LFUQ-WH1, we used a similar approach to Sayer et al. (2010), correlating the core with LFUQ-M1 according to specific independent paleo-markers: changes in LOI profiles, changes in lithostratigraphy, changes in geochemical proxies, and changes in pollen abundance of Pinus spp. which was introduced in the lake catchment during the early 1900s (Vidal and Perez, 2007).

We analyzed 25 (LFUQM1) and 19 (LFUQB1) sediment samples (1-cm thick) for plant macro remains at 2-cm intervals for the top 20 and at 5 cm intervals lower down each core, respectively. Between 20 and 40 cm³ of wet sediment was disaggregated in 10% potassium hydroxide (KOH) before sieving. Macrophyte remains were retrieved from the residues of sieved core material (using mesh sizes of 355 and 125 μm) following standard methods (Birks, 2001). Plant remains included seeds and fruits, leaf-spines, leaf fragments (including water lily leaf sclereid cells), charophyte oospores and characteristic Eichhornia root-caps. Plant remains data were standardized as the number of fossils per 100 cm³ and identified by comparison with reference material and by using relevant taxonomic keys (e.g., Birks, 2001; https://tmi.laccore.umn.edu/is). Due to poor macrofossil preservation, we estimated temporal variations of Egeria via pollen analysis using the method for determining absolute pollen concentration (Moore et al., 1991). Fifteen sediment samples were analyzed for pollen analysis for LFUQ-M1 core and 12 samples for LFUQ-WH1.

We assessed relationships between current macrophyte species occurrences and variation in physical and chemical parameters across transects via multivariate regression tree (MRT) analysis (De'ath and Fabricius, 2000). Prior to analysis, we excluded variables with collinearity (i.e., conductivity and TDS; see Table S1.1 in Appendix S1), defined as those with a variance inflation factor >8 (Zuur et al., 2007). To match the number of environmental data points at each transect, we split the macrophyte data from each transect into 3 groups according to proximity to the beginning, middle or end of each transect. The occurrences of each macrophyte species within each group transect were then transformed as frequencies of occurrence (i.e., [number of observations/total number of sampling points per group transect] X 100), and environmental data were log(x+1) transformed prior to analysis (Legendre and Gallagher, 2001). Cross-validation was applied to determine the optimal number of splits in the tree (De'ath and Fabricius, 2000) and the best predictive tree with one standard error was selected (De'ath and Fabricius, 2000). Given that the MRT algorithm uses random numbers as a starting point, we ran several (n = 15) MRT analyses using the set.seed algorithm in R (R Development Core Team, 2016) and compared the multiple runs for the same model (De'ath and Fabricius, 2000). The mean value of the best predictive trees was then used to describe the variation in the macrophyte and environmental variables. MRT analysis was carried out using the MVPART library (version 1.4; De'ath and Fabricius, 2000) in R (R Development Core Team, 2016). Principal component analysis (PCA) was then used to visualize macrophyte species distributions according to the MRT selected groups.

To reduce the variation between differentially produced plant macrofossil structures (e.g., root-caps, leaves and seeds), we log(x+1) transformed the plant macrofossil abundance (Legendre and Gallagher, 2001). In order to assess major temporal zones of community change, we combined the data sets from the two cores into a single species matrix and applied a univariate regression tree analysis using core depth as the explanatory variable. To explore any potential variability between different tree models, we ran 10 different analyses using the set.seed algorithm in R (R Development Core Team, 2016). The mean value of the best predictive trees was then used to describe the variation in the plant macrofossil data. A PCA was used for visualization of major compositional temporal changes in macrophyte assemblages across the cores. Plant macrofossil stratigraphy for each core was achieved using the “Rioja” Package in R (Juggins, 2009; R Development Core Team, 2016).

To assess long-term dynamics of invasive macrophyte species spread, we focused on the variation of Eichhornia fossil abundances (longer history of occurrences and good fossil preservation). We modeled the variation of Eichhornia log-abundances through time as a two distinct spread processes using piecewise regression analysis (Aikio et al., 2010). This analysis models early stages of spread (i.e., stable periods of low abundances; lag-phase) as a linear process while the subsequent exponential growth phase is modeled as a non-linear component (Aikio et al., 2010). The significance of a potential breaking point between the two-modeled slopes was tested using a Davies test (Davies, 1987). Analyses were run using the “fragmented” package in R (R Development Core Team, 2016; Muggeo and Muggeo, 2017).

A total of 14 macrophyte species and filamentous algae (undifferentiated) were recorded across the four studied transects including Azolla filiculoides, E. crassipes, E. densa, Limnobium laevigatum Humb. & Bonpl. ex Willd. Heine, Ludwigia peploides Kunth P.H.Raven, Myriophyllum aquaticum (Vell.) Verdc., Hydrocotyle ranunculoides L. f. (VVV), Potamogeton pusillus L., Utricularia gibba L., Typha angustifolia Pers., and Schoenoplectus californicus (C.A. Mey.) Steud (Table 3). Plantless areas were a common feature (mean = 38.9% ± 31.9) across the open water transect (OW1) and, when present, E. densa dominated the assemblages (mean = 51 ± 34.8%). Other species such as E. crassipes (mean = 20 ± 26.1%), P. pusillus L. (mean = 6.7 ± 15.2%), and A. filiculoides (mean= 5.6 ± 8.1%) occurred irregularly. In the Eichhornia type transect (WH1) we recorded 14 species including E. crassipes, L. laevigatum, L. peploides, M. aquaticum, H. ranunculoides, and U. gibba. Emergent species like Bidens laevis (L.) Britton Sterns Poggenb., T. angustifolia, and S. californicus, also occurred in a few sample locations. In ED1 (E. densa type) transect Egeria dominated and a few species such as A. filiculoides, L. laevigatum and U. gibba occurred sporadically (Table 2). For the transect M1 (mixed assemblage of species) we observed 10 macrophyte species with E. densa (100%), A. filiculoides (mean = 97.8 ± 8.6%) and Eichhornia (mean = 78.9 ± 18.3%) being the most common ones.

A summary of measured environmental variables is presented in Table 4. Transects OW1 (137 ± 36.5 cm), M1 (114.8 ± 48 cm) and WH1 (103 ± 41.4 cm) were generally deeper than the ED1 (86 ± 31.1 cm) transect. In a descending order, DO surface values varied from 12.3 ± 2.9 ppm at OW1, 11.8 ± 2.1 ppm at ED1, 11.4 ± 2.8 ppm at M1 and 6.7 ± 2.6 ppm at WH1. Similarly, DO in the water column had the lowest values at the WH1 (3.9 ± 2.6 ppm) transect. Highest secchi depths were recorded at the E1 (55 ± 17.9 cm) and M1 (50.8 ± 11.7 cm) transects, while lowest values (19.3 ± 6.1 cm) were observed at the WH1 transect. Salinity (0.25 ± 0.05 ppt), conductivity (357 ± 99.8 μS/cm) and TDS (242 ± 62.6 ppm) were also greater at the ED1 transect than at the other three transects.

The MRT analysis identified surface DO, secchi disk depth and salinity as the most important environmental variables explaining the current macrophyte species distributions across the transects (Figure 2A). Nine out of the 15 analyzed trees split the macrophyte data into five main groups explaining 67.3 ± 3.4% of the plant variation (Figure 2A; Table S2.1 in Appendix S2). PCA on the macrophyte frequency of occurrences showed that the MRT tree splits generally concur with different macrophyte cover types of each sampling transect (Figure 2B). The tree was initially split in two main branches, where the majority of macrophyte data points of transect WH1 clustered in an independent group (lower right hand side in Figure 2B) characterized by DO surface values <7.6 ppm. Three sampling points characterized by Eichhornia dominance clustered, as an independent group with DO values above 7.6 ppm. Macrophyte data points from transect OW1, ED1, and M1 grouped together on the left branch of the tree according to DO surface values ≥7.7 ppm (Figure 2A). The MRT separated most macrophyte data associated with the OW1 transect according to DO water column values ≥ 7.73 ppm, whereas higher salinity for ED1 (≥0.15) and higher secchi depths (≥30 cm) split most ED1 data from M1 data. The PCA analysis showed that macrophyte associations of these three groups overlap in some instances (Figure 2B).

Our age-depth calibration of the core LFUQ-M1 resulted in a chronological model spanning the last c.600 years, with the last two centuries in the top 20 cm (Figure 3; Table 2). Due to the non-monotonic variation of unsupported ²¹⁰Pb in the core, radiometric chronologies were calculated using the CRS dating model (See Appendix S4 for details). The simple CRS model placed the 1963 depth above 17 cm, likely due to low unsupported ²¹⁰Pb activities and to earlier reductions in unsupported ²¹⁰Pb. Chronologies and sedimentation rates were therefore corrected under the CRS model by assigning a date of 1963 to sediments at 19 cm in accordance to the ¹³⁷Cs record. The ¹⁴C AMS analysis dated the macrofossil samples at 30, 75, and 98 cm to 1676 ± 30, 1506 ± 40, and 1370 ± 86, respectively (Table 2). Sedimentation rates in our age-depth calibration showed a U-shaped pattern being relatively high (1 cm/yr.) between 98 and 76 cm, dropping to a mean of 0.0065 cm/yr between 75 and 20 cm and gradually increasing above 20 cm to reach a mean value of 2.5 cm/yr in the most recent samples (6–1 cm) (Figure 4).

Temporal variations in lithostratigraphy, LOI and XRF index ratios were highly consistent across the two cores, showing three zones of change that allowed us to correlate the dated core LFUQ-M1 with the LFUQ-WH1 core with good confidence (Figure 4). The three zones correspond to: (i) a historical phase with distinctive dark sediments, high organic matter, high Ca/Ti, (Fe/Mn)/Ti and Rb/K ratios, and low K/Ti ratio (100–75 cm in LFUQ-M1; 50–41 cm in LFUQ-WH1); (ii) a transitional phase with distinctive pale gray sediments, low organic material, low Ca/Ti, (Fe/Mn)/Ti and Rb/K ratios, and hig K/Ti ratio (74–20 cm in LFUQ-M1; 38–13 cm in LFUQ-WH1); (iii) a recent phase with distinctive dark sediments, high organic matter, Ca/Ti, (Fe/Mn)/Ti and Rb/K ratios and poor K/Ti ratio (19–0 cm in LFUQ-M1; 13–0 cm in LFUQ-WH1 cm). Similarly, the independent pollen record of Pinus spp. was highly consistent across the two cores showing a reduction at 9–10 cm, an increase at 7.5 and a reduction at 2 cm (see Figure S5.1 in Appendix S5).

The plant macrofossil analysis resulted in 14 macrophyte taxa occurring throughout the cores (Figures 5, 6). Plant assemblage representation across the two cores was highly congruent resulting in nine out of 10 consistent regression trees explaining the macrophyte variation of 53% ± 4.1 of (Figure 7; Table S3.1 in Appendix S3). Three major temporal zones of change were identified: Zone 1– CE c. pre-1506 (100–75 cm in LFUQ-M1; 50–41 cm in LFUQ-WH1); Zone 2– 1506-1904 CE (74–20 cm in LFUQ-M1; 38-13 cm in LFUQ-WH1); and Zone 3– 1905 CE-present (19-0 cm in LFUQ-M1; 13-0 cm in LFUQ-WH1 cm). In Zone 1 charophytes (Chara sp. and Nitella sp.), bryophytes (agg.), Ceratophyllum demersum L., and P. pusillus, dominated the assemblages. Leaf-sclereid cells (probably Nymphaea sp.), Potamogeton illinoensis Morong., the previously unreported Najas guadalupensis (Spreng.) Magnus, A. filiculoides and E. crassipes also occurred in abundant numbers. In Zone 2 charophytes, bryophytes and C. demersum disappeared from the fossil record, while Nymphaea spp., P. illinoensis, and N. guadalupensis dominated until the end of the 1800s when abundances gradually decreased. Eichhornia and A. filiculoides occurred in relatively low abundances, presenting a slight increase from 1640. At zone 3 (1905-present), Nymphaea sp., disappeared while N. guadalupensis and P. illinoensis further declined. At the beginning of this zone E. densa colonized the lake whilst E. crassipes and A. filiculoides gradually increased. The species L. valdiviana and M. aquaticum also become apparent during the late periods of this zone and since the mid-1990s the lake was typified by marked expansions in Eichhornia, A. filiculoides, Egeria, M. aquaticum and L. peploides. At this point P. illinoensis disappeared from the fossil record while N. guadalupensis only occurred in low numbers throughout the LFUQ-WH1 core (Figure 6).

The palaeo-data suggest that both Eichhornia and Egeria have been present in Fúquene Lake well before the proposed date of introduction in the 1950s and 1990, respectively (Figures 5, 6). Our data suggest that Egeria colonized the lake around the early 1900s, while Eichhornia appeared at 1500 CE. Piecewise regression analysis on Eichhornia temporal variation showed two phases of spread explaining 86.4% of the adjusted variation (_adjR²) with a highly significant (P< 0.0001) breaking point at 9.2 ± 0.6 cm, which corresponds to 1995–2000 (Figure 8). Accordingly, the Davies test revealed a highly significant (P < 0.0001) change in slopes at 11.5 cm (1995) characterized by large increases in fossil abundances that almost doubled the historical values (Figure 8). These patterns are highly consistent with the historical records of Eichhornia spread in the lake (Figure 1; Table 1).

Our combined biological and geochemical palaeolimnological approach showed that over the last c.600 years macrophyte community dynamics in Fúquene Lake have gradually changed in response to long-term natural and human induced variations in lake catchment dynamics, hydrological modifications (lake size), reductions in water quality and the invasion of macrophyte species. For instance, the high Ca/Ti ratio suggests that pre-1500 climatic conditions were dry and the mean lake water balance was negative (Haberzettl et al., 2009). Shallower lake conditions would stimulate lake productivity (high LOI) and associated reductive conditions in the lake sediments (high [Fe/Mn]/Ti ratio). Chemical weathering (high Rb/K) from exposed catchment soils would have been more intense (Davies et al., 2015). The occurrence of sensitive macrophyte species such as charophytes (Chara sp. and Nitella sp.), bryophytes (agg.), and P. illinoensis further suggests oligo-mesotrophic conditions at the time (Kolada et al., 2014).

Variations in the geochemical composition after 1500 CE suggest an environmental shift toward a wetter climate and lake expansion (low Ca/Ti; Haberzettl et al., 2009). This climatic pattern was described by Bird et al. (2018) in a nearby neotropical Andean lake (Ubaque Lake, Colombia). As indicated by the high K/Ti ratio (Figure 4), physical erosion may also have increased during this time period in response to larger river inputs and human activities in the catchment area (Arnaud et al., 2012). Macrophyte compositional shifts at this time suggest incipient signs of nutrient-enrichment and/or reductions in water transparency with the disappearance of charophytes and bryophytes from the fossil record. These two macrophytes are commonly associated with low-nutrient, clear waters and decrease as more nutrient-enriched and turbid waters develop (Poikane et al., 2014; Ecke, 2018). Although archaeological and paleoecological records from the lake catchment area have shown that humans occupied and modified the landscape for agriculture well before 1500 CE (Van Geel and Van der Hammen, 1973; Langebaek, 1995) it is likely that this change of lake water quality is linked to major shifts in agriculture and increases in population densities after the Spanish conquest at 1500 CE (Langebaek, 1995). The inferred wetter conditions during this colonial time period, together with technological changes in agricultural (large numbers of grazing livestock) and mining practices initiated by European arrival (Etter et al., 2008) would have resulted in greater sediment and nutrient runoff into the lake. Compositional shifts in lake biota associated with the onset of European settlement have been similarly reported across subtropical lakes (e.g., Bickford et al., 2008; Reid, 2008; Reid M. A. et al., 2018).

Abundances of some existing macrophyte species (Nymphaea sp., N. guadalupensis, P. illinoensis, E. crassipes, and A. filiculoides) initially increased after 1500 CE with the land use transformations and wetter conditions (Figures 5, 6). However, N. guadalupensis, P. illinoensis and Nymphaea sp. started to decline around the late 1800s. By the mid 1800's CE, agricultural areas had doubled in size (Etter et al., 2008), likely resulting in greater nutrient load into the lake. Further reductions of N. guadalupensis, P. illinoensis, and Nymphaea sp. after the early 1900s CE coincides with the lake drainage schemes (as supported by the increases in Ca/Ti and [Fe/Mn]/Ti ratios) which aimed to further expand grazing land (Figure 1; Franco et al., 2011). Such hydrological alteration would have resulted initially in turbid waters coupled with heavier nutrient loads (Table 1), which are unfavorable for N. guadalupensis, P. illinoensis, and Nymphaea sp. (Wentz and Stuckey, 1971; Skubinna et al., 1995; Egertson et al., 2004). In particular, N. guadalupensis, and P. illinoensis have been reported to expand under mesotrophic conditions (Wentz and Stuckey, 1971; Egertson et al., 2004), and to be extirpated as turbidity and/or nutrient-enrichment progresses up to eutrophic levels (Skubinna et al., 1995; Brush and Hilgartner, 2000; Egertson et al., 2004).

The marked expansion of Eichhornia, Azolla and other associated taxa were most likely the result of significant increases in nutrients (0.41 mg/L by 2010) and further shifts in lake hydrology (Table 1; Figure 4). The eutrophic conditions coupled with the implementation of the surrounding water canal would have created ideal conditions for floating plants to thrive, in particular Eichhornia and Azolla (Scheffer et al., 2003; Villamagna and Murphy, 2010). The high capacity of E. densa to obtain nutrients from both the sediments and the water column, and to grow under low light conditions (Barko and Smart, 1981; Yarrow et al., 2009) may have facilitated its spread. Population growth, along with the deliberate introduction of new Egeria plant propagules in 1990, would have helped to out-compete the already declining populations of N. guadalupensis and P. illinoensis.

Contrary to previous suppositions, our data from two cores showed that E. densa and E. crassipes have been present in Fúquene Lake well before the suggested dates of introduction (Figures 5, 6). According to herbarium collections, E. densa has been present in the Colombian Eastern Andean Cordillera since at least the mid-1950s (Schmidt-Mumm, 1996) and Fúquene Lake technical reports suggested its presence since 1986 (JICA-Agencia de Cooperación Internacional Del Japón, 1999). Although the origin of E. crassipes is probably Brazilian Amazonia (Barrett and Forno, 1982), the species spread north, and had colonized many countries of Central America and the Caribbean by the end of the nineteenth century (Barrett and Forno, 1982). It is therefore probable that these two species occurred in Fúquene well before current assumptions but given its low abundances they might have been overlooked (Pyšek and Hulme, 2005).

Our data showed that both Eichhornia and Egeria demonstrated classic invasion behavior (sensu Blackburn et al., 2011): they had already crossed an environmental barrier to spread, they showed a distinct lag-phase before the 1990s CE (low abundances) and grew exponentially to high levels after the 1980s (Figure 8). This bi-phasic temporal behavior shows that during moderate nutrient conditions (pre-1930s), both invasive species coexisted with native aquatic plants, but when conditions became eutrophic they thrived. Eichhornia occurs as a free-floating or, to a lesser extent, as a submerged or emergent plant (Gopal, 1987). Thus, it can endure extreme habitat disturbances such as those observed in Fúquene Lake (Gopal, 1987). The lake surface area reduction during the 1930s may have resulted in large, exposed, littoral areas and shallower waters (as suggested by the high Ca/Ti and Rb/K ratios) with stable daily and seasonal average water temperatures, all of which have shown to stimulate seed germination and plant growth (Wilson et al., 2005; Albano Pérez et al., 2011). Ultimately, the more steady hydrological conditions and the heavy nutrient load derived from the surrounding canal constructed the mid 1980s would have stimulated Eichhornia proliferation to the current dominant levels (Figure 1; Reddy and DeBusk, 1984; Villamagna and Murphy, 2010). Thus the loss of resident macrophyte species may not be fully accounted for by the spread of invasive Eichhornia and Egeria but through habitat degradation that decreased the competitive capacities of resident species such as Nymphaea sp., N. guadalupensis and P. illinoensis whilst simultaneously assisting the spread of two invasive macrophytes (MacDougall and Turkington, 2005). The gradual increases in (Fe/Mn)/Ti since the mid-1900s along with the contemporary monitoring data indicate that, once fully dominant, these two invasive species along with Azolla may have been exerting secondary effects on the remaining biota via reductions in DO, intensification of light limitation in the case of Eichhornia and Azolla (Table 2; Albright et al., 2004; Villamagna and Murphy, 2010) and via competitive exclusion in the case of Egeria (Yarrow et al., 2009).

We are aware that the use of paleolimnological data to infer past communities has limitations. Due to preservation and a strong likelihood of missing rare or distantly-located macrophyte taxa (Zhao et al., 2006), not all species present historically will leave remains in sediment cores. Nevertheless a substantial portion of current-day macrophyte species (in this case 60%) was represented in the surface sediments of the two sediment cores (in agreement with Davidson et al., 2005; Salgado et al., 2010). Uncertainties in the chronology model, particularly for dates below the ²¹⁰Pb record, might have introduced some discrepancies regarding the exact dates of invasive species colonization (Figure 3; Table 2). Nonetheless, our data showed that E. densa and E. crassipes were present in Fúquene Lake well before the suggested dates of introduction when dating errors are accounted for.

Our study suggests that once invasive species dominate, they exert positive (e.g., enhancing lake macrophyte biodiversity) and negative (e.g., reductions in DO) secondary effects on the remaining resident submerged species, and influence variability in environmental parameters across the lake. For instance, at the lake scale, the clear spatial separation between areas dominated by Egeria and Eichhornia and their unique associations with other macrophyte species and limnological variables suggest that these species are currently creating significant lake habitat heterogeneity. Such spatial patterns seem to be helping to maintain a relatively diverse macrophyte community, at least in the short-term, even under the current highly eutrophic conditions (Figure 2). These patterns in macrophyte community, environmental heterogeneity and species richness may help to increase lake resistance to nutrient enrichment (Olden et al., 2004; Salgado et al., 2018b).

Egeria areas had low macrophyte diversity and were highly homogenous in macrophyte composition. These patterns confirm previous studies on Egeria that have shown that it outcompetes other submerged macrophytes and so, dominates under stressful conditions (Roberts et al., 1999; Yarrow et al., 2009). Competitive adaptations include fast-growth rates, aggressive vegetative propagation, and a capacity to grow under low light conditions (Barko and Smart, 1981; Yarrow et al., 2009). Egeria is apparently a suitable plant for phytoremediation because it increases DO under the current eutrophic conditions (Table 3). As a submerged species this plant releases O₂ directly into the water column, thus helping increasing oxygen levels (Yarrow et al., 2009). Higher secchi depths across dense stands of Egeria relative to the other transects, also supports the idea that stands of this species help to decrease water movement, reduces sediment resuspension and improves water transparency (Yarrow et al., 2009). Egeria forms denser and more spatially extensive plant stands than the struggling native submerged N. guadalupensis, and P. illinoensis, (Yarrow et al., 2009). Therefore, its dominance may be providing greater habitat opportunities to associated faunal organisms (e.g., fish and invertebrates) that depend on macrophytes for food or as a substrate.

As documented elsewhere (Villamagna and Murphy, 2010), we found low levels of DO and reduced water clarity across Eichhornia areas. Unlike Egeria, the dense floating mats of Eichhornia releases oxygen into the atmosphere while preventing light reaching the water column (Villamagna and Murphy, 2010). Thus, anoxic conditions and light limitation beneath the floating mats constrain the survival of phytoplankton and submerged macrophytes such as N. guadalupensis and P. illinoensis which are poorly adapted to low light conditions (Davis et al., 1985; Brush and Hilgartner, 2000; Egertson et al., 2004). Eichhornia areas appear to be macrophyte rich (Figure 2). The floating mats create ideal island-like conditions for other emergent plant species, in this case B. laevis, L. peploides, S. californicus, and T. angustifolia, to grow (Gopal, 1987); These free-floating mats also create a series of temporary open water inlets and semi-enclosed basins with low wave exposure offering temporal refugia for other floating-plants such as A. filiculoides, H. ranunculoides, L. laevigatum, and also submerged species like M. aquaticum, E. densa, P. pusillus and, U. gibba (Figure 2).

Studies on the effects of eutrophication on neotropical/subtropical lake macrophytes suggest that nutrient enrichment reduces the resilience of the system against a shift to a state dominated by floating plants (Scheffer et al., 2003). In accordance, Eichhornia, L. laevigatum, L. peploides, and Azolla currently cover large areas of Lake Fúquene. Our data do not suggest a sudden tipping point toward a shift in floating plant dominance. Instead, we see gradual transitions in composition through time (Figures 5, 6) where, despite large nutrient load and associated increases in turbidity, the extensive floating-plant mats still co-exist with dense stands of submerged plants (Figure 2). Co-existence of floating and submerged plant species is largely due to the recent spread of Egeria, highlighting that alternative regime shifts to a floating-plant dominated state may not fully depend on nutrient availability but also on the competitive ability of submerged species to endure the light-limited conditions imposed by eutrophication. The transition between submerged and floating macrophytes over time provides support for the growing recognition that in some natural shallow lakes, macrophyte communities may not change abruptly between ecological regime shifts (Scheffer et al., 2003) but gradually over long-term scales in response to eutrophication (Sayer et al., 2010).

To date, there are few studies addressing the full extent of abiotic and biotic changes that may occur in response to the spread of invasive macrophyte species in neotropical lakes (Villamagna and Murphy, 2010; Rejmánková et al., 2018). Even less is known of the ecological conditions prior to invasion and the long-term dynamics of spread linked to other co-occurring stressors such as eutrophication (Strayer et al., 2006). Our results suggest that equating macrophyte biodiversity loss with the spread of invasive macrophytes in degraded lakes could be misleading. Instead, we found that habitat degradation is a major driving force of biodiversity loss and invasive species spread. However, once invasive macrophyte species dominate, they may exert both synergistic and antagonistic secondary effects on plant biodiversity by modifying habitat structure (e.g., shading and anoxic conditions), by excluding remaining resident species through direct resource competition (e.g., Egeria) and/or by exerting habitat heterogeneity across the lake.

Our data show that Fúquene Lake is in a delicate hybrid ecosystem state (sensu Hobbs et al., 2009), where, despite the depletion of most native submerged species, the lake still retains some historically characteristic species like E. crassipes, A. filiculoides, and P. pusillus (Figures 5, 6). The current dominance and function of Eichhornia and Azolla lies outside its historic range of variability and is working synergistically with eutrophication. Furthermore, Eichhornia, Azolla and Egeria are extremely difficult to eradicate once they dominate, and any further nutrient inputs will encourage their spread across the lake. Thus, restoring the lake to a state closer to its historical conditions is likely impossible and/or economically impractical (Hobbs et al., 2009). We therefore suggest focusing management objectives on reducing nutrient inputs and further controlling the spread of the two invasives along with Azolla to enhance open water areas that may potentially allow the currently threatened historical submerged macrophytes N. guadalupensis and P. illinoensis to re-grow. Addressing the way the surrounding canal is currently managed and restoring the lateral and longitudinal hydrological connectivity are likely key to accomplishing these conservation goals.

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.