Hydrological alteration, eutrophication and macrophyte invasion generate novel ecological states and biotic homogenisation in shallow lakes. As these stressors increasingly co-occur and unfold over decades to centuries, disentangling their combined effects on invertebrate assemblages remains however, challenging. To assess their long-term (decades-centuries) interactions on driving aquatic invertebrates, we analysed three ²¹⁰Pb-dated sediment cores spanning contrasting dominance of the invasives Egeria densa (submerged), Pontederia crassipes (floating) and Azolla filiculoides (floating) in Fúquene, a largely-drained, eutrophic Andean shallow lake (Colombia). Eighteen invertebrate taxa were recorded with multivariate analyses explaining 51% of assemblage variation and identified a pre-1800 mesotrophic phase; a drainage phase (1800–late 1980s) marked by shifts towards eutrophic taxa; and a post-perimeter canal phase (post-1990) characterised by invasive expansion and distinct macrophyte-associated invertebrate assemblages. Despite restructuring, invertebrate assemblages did not homogenise but became more heterogeneous among cores in Zone 3. Variation partitioning revealed strong coupled hydro-chemical–macrophyte effects at the assemblage and taxon level explaining 13–56% of the variation, with invasive growth forms uniquely accounting for 5–40%. E. densa stands supported aquatic insects, whereas floating species favoured Oribatidae, bryozoans and semi-terrestrial chironomids. Invasion–eutrophication-driven homogenisation is therefore not inevitable but contingent on hydrology, environmental variation and growth-form diversity.

1.   Material and methods

(a) Study area

Lake Fúquene is a shallow, tropical, upland lake (2,560 m elevation) in the eastern Cordillera of the Colombian Andes (5°27'55''N, 75°46'19''W) (figure 1). With a current surface area of 30 km² and mean depth of 2 m (maximum 6 m), the lake is eutrophic (total phosphorus [TP] > 0.40 mg/L; total nitrogen [TN]: 1.8 mg/L), with TP and TN loading predominantly from agriculture (60%) and domestic sources (~20%) [55,56]. TP and TN concentrations in surface sediments are also elevated (TN: 2,095 mg/kg; TP: 815.08 mg/kg) [56]. Since the 1990s, the lake has been dominated by large stands of P. crassipes, A. filiculoides, and E. densa, creating distinct limnological conditions (figure 1). Areas dominated by P. crassipes and A. filiculoides support higher macrophyte diversity (10-14 species) but lower dissolved oxygen (DO < 4 mg/L) and water transparency (< 50 cm). In contrast, E. densa-dominated areas present lower macrophyte richness (5 species) but higher DO (> 6 mg/L), transparency (> 55 cm), and conductivity (> 300 µS/cm) [28].

(b) Palaeolimnological data

To quantify temporal changes in invertebrate assemblages and the spread of three invasive macrophytes, we retrieved three short sediment cores (<1 m) from: (i) an area dominated by P. crassipes (core code: LFUQ-WH1), (ii) an area dominated by A. filiculoides with mixed patches of P. crassipes and E. densa (LFUQ-M1), and (iii) an area dominated by E. densa (LFUQ-ED1) (figure 1). The cores were collected in 2016 during an unusually dry El Niño–ENSO event, which in Fúquene Lake has been shown to drastically reduce water levels, enhancing both floating and submerged macrophyte cover [28;56] (figure 1). We used a large-bore corer (10 cm diameter, Aquatic Instruments Inc.) to secure substantial sediment volume while preserving undisturbed surface sediments [57]. To minimise potential variations from differences in proximity to shore and water depth, all cores were collected from the western side of the lake at similar depths: LFUQ-M1 (5°27'13.86"N, 73°44'59.33"W; 100 cm long) and LFUQ-WH1 (5°27'45.35"N, 73°45'48.61"W; 50 cm long) at 100 cm depth, and LFUQ-ED1 (5°27'45.86"N, 73°45'58.98"W; 50 cm long) at 90 cm depth. All cores were extruded in the field at 1 cm intervals, with lithostratigraphic changes recorded prior to extrusion.

(i) Palaeo-analysis

We analysed 56 sediment samples (1 cm thick) across the three cores: 25 from LFUQ-M1, 16 from LFUQ-WH1, and 15 from LFUQ-ED1, obtained at 2 cm intervals from the top 20 cm and 5 cm intervals thereafter. Invertebrate and plant macrofossil relative abundances analyses followed Salgado et al. [21]. For each sample, 20-40 cm³ of wet sediment was disaggregated in 10% KOH and sieved through 355 µm and 125 µm meshes. To process invertebrates and macrofossils simultaneously, we modified standard chironomid methods by using 125 µm sieves instead of conventional 90 µm [21], potentially missing smaller taxa. Invertebrate remains included Oribatidae (body fragments), Bryozoa (Plumatella spp., statoblasts), Daphniidae (Daphnia spp., Ceriodaphnia spp., Simocephalus spp., ephippia), Decapoda (body fragments), aquatic insects (Diptera, Corixidae, Odonata body fragments), and Chironomidae (head capsules). Chironomid head capsules were mounted and identified using standard methods [58] to genus or type level using reference collections and keys [58-60]. All relative abundances were standardised as fossils per 100 cm³. Plant macrofossil data (A. filiculoides, P. crassipes) and E. densa pollen from LFUQ-WH1 and LFUQ-M1 were obtained from Salgado et al. [28]. For LFUQ-ED1, lacking pollen data, temporal variations were established through correlation with the other cores and historical records [28].

            As eutrophication proxies, we used loss-on-ignition (LOI) [61] to determine changes in organic carbon fraction in core sediments, along with TN data for the LFUQ-M1 core. TN was measured by combustion of dry sediment on a Costech Elemental Analyzer coupled to a VG TripleTrap and Optima mass spectrometer at the National Environmental Isotope Facility, British Geological Survey, UK. Although sedimentary TP data were unavailable, TN use as a eutrophication proxy is well-supported by historical and surface sediment records showing TN and TP increased in parallel with eutrophication [56]. Moreover, palaeolimnological studies in other Colombian shallow lakes have demonstrated that TP strongly correlates with TN and LOI [62].

            As a proxy for hydrological change (drainage and perimeter canal implementation), we used the elemental ratio of zirconium to iron (Zr/Fe) [62]. This ratio has been shown to be a reliable proxy for hydrological energy variation in Colombia [62,63], with coarse-grained Zr-rich sediments deposited during periods of high energy from floods, while finer Fe-rich sediments prevail during dry periods or, in our case, following lake drainage and canal implementation. These geochemical elements were obtained from Salgado [28] and measured on 1 cm-thick discrete dried, homogenised core samples using a portable XRF analyser spectrometer (XMET 7500).

(ii) Dating and age-depth models

The age model for LFUQ-M1 was derived from Salgado et al. [28] using radiometric measurements of ²¹⁰Pb, ²²⁶Ra, ²¹⁴Pb and ¹³⁷Cs, applying a Constant Rate of Supply (CRS) model [64], and supplemented by ¹⁴C AMS analysis of plant macrofossils (electronic supplementary material, Table S1). In the absence of direct radiometric dating for LFUQ-E1 and LFUQ-WH1, chronologies were established by correlating geochemical proxy profiles (LOI, Zr/Fe, K and Ti) with LFUQ-M1 using multivariate Dynamic Time Warping (DTW) implemented in the dtw package in R [65]. Proxy series were interpolated to a common depth grid, square-root transformed where appropriate, and standardised prior to alignment. DTW was applied using Euclidean distance with a symmetric step pattern and a Sakoe–Chiba window constraint to preserve stratigraphic coherence.

            To account for uncertainty, multiple DTW models were generated using different proxy combinations. Age estimates were summarised as minimum, maximum, median, and mean values, defining uncertainty envelopes. Final age–depth models were selected based on agreement among models, stratigraphic consistency, and minimisation of age reversals. Mean ages were adopted, with minor adjustments to maintain monotonic depth–age relationships and consistency with historical records. Although an age model for LFUQ-WH1 was previously derived through visual multi-proxy stratigraphic correlation by Salgado et al. [28], we recalibrated this chronology using the same DTW-based framework applied to LFUQ-E1 to ensure methodological consistency across cores.

(c) Data analysis

We analysed the invertebrate palaeo-record in relation to three main historical environmental stressors: (i) hydrological change (Zr/Fe), (ii) eutrophication (TN), and (iii) macrophyte variation. Note that LOI was excluded from the data analysis due to its high collinearity with TN (r = 0.72) and Zr/Fe (r=-0.82) (electronic supplementary material, figure S1) but kept for discussion.

            To assess whether macrophyte variation influenced invertebrates via overall assemblage shifts (i.e., changes in composition and abundance of native and non-native species) or through changes in specific invasive taxa, we partitioned the plant macrofossil data into three predictors: (1) overall macrophyte assemblage structure, (2) floating invasives (P. crassipes and A. filiculoides), and (3) submerged invasives (E. densa). Overall assemblage structure (1) shifts at each core were summarised using principal curve analysis implemented in the analogue package in R [64].

      All environmental and ecological variables were square-root transformed prior to analysis to reduce extreme variation (~300-fold) and down-weight dominant values [65]. Stratigraphic plots were produced using the tidypaleo [67] and ggplot2 [68] packages in R and major phases of ecological change were identified using stratigraphically constrained hierarchical clustering (CONISS) applied to square-root transformed invertebrate assemblage data in the rioja package [69]. The square-root transformation was used to down weight dominant taxa and stabilise variance prior to clustering. As CONISS operates on Euclidean distances, clustering was performed on the transformed data under a stratigraphic constraint that preserves sample order. The optimal number of zones was evaluated using a broken-stick stopping rule applied to dendrogram structure. Given the progressive nature of ecological change, zone boundaries were interpreted as transitional and defined based on dendrogram structure rather than fixed chronological thresholds.

(i) Gradients of ecosystem change

We applied Multiple Factor Analysis (MFA) to identify and simultaneously summarise multivariate temporal gradients in invertebrate assemblage composition, eutrophication (TN), hydrological change (Zr/Fe) and macrophyte variation (overall assemblage structure and invasive floating and submerged) [72]. The invertebrate taxa were grouped into Oribatidae, Bryozoa (Plumatella spp.), Daphniidae (Daphnia spp., Ceriodaphnia spp. and Simocephalus spp.), Decapoda, aquatic insects (Diptera, Coxidae and Odonata), and Chironomidae (morphotypes). The square root-transformed data were standardised during analysis to create comparable units, and the MFA was conducted MFA using the FactoMineR [73].

            We then assessed whether invertebrate assemblages (composition and relative abundances) became more homogeneous/heterogeneous over time in response to major stressors using tests of homogeneity of multivariate dispersion (HMD) [74] implemented with betadispers in vegan [70]. Groups were defined according to the principal stratigraphic zones identified by CONISS, representing the main phases of assemblage change in the palaeo-record. Euclidean distances were calculated from individual sample scores (first five axes, >70% of total variance) derived from a second MFA including only invertebrate data. Mean distances to group medians were compared using Tukey’s Honest Significant Difference test (p < 0.05), where higher values indicate greater assemblage heterogeneity and lower values indicate homogeneity. Because sample numbers differed among CONISS-defined groups, distances were adjusted using the n/(n−1) correction to account for unequal sample sizes.

(ii) Independent vs. joint effects of environmental drivers

To quantify how environmental drivers structured invertebrate communities, we applied partial Redundancy Analysis (pRDA) at both assemblage and taxon levels in vegan [70]. Environmental predictors were grouped into a hydro-chemical component (TN, Zr/Fe) and a macrophyte component (total assemblage, floating invasive, submerged invasive). Assemblage-level analyses were used to detect lake-scale patterns in community composition, whereas taxon-level analyses targeted between-core variability associated with differences in macrophyte growth forms.

            At the assemblage level, forward selection was applied to square-root-transformed community data to identify significant predictors, retaining variables that improved model fit (ordistep, 999 permutations, p < 0.05). The selected predictors were then used in pRDA, and variation partitioning (varpart) was applied to decompose explained variance into unique and shared components. The significance of unique fractions was tested using partial RDA (999 permutations, p < 0.05). At the taxon level, forward selection was similarly performed prior to pRDA using RDA from a null model (taxon ~ 1), and retained predictors were incorporated into pRDA models to assess species-specific responses following the same approach as for the assemblage-level analysis.