Biodiversity Dynamics in a Ramsar Wetland: Assessing How Climate and Hydrology Shape Distribution of Dominant Native and Alien Macrophytes

1. Introduction

Wetlands, and particularly coastal wetlands, fulfil numerous ecosystem functions and deliver essential ecosystem services to human society across the world. These include the preservation of native biodiversity, regulation of biogeochemical processes and water cycles, carbon sequestration, provision of habitat for commercial fisheries, storm surge protection, and opportunities for cultural, recreational, and educational activities [1,2,3,4,5,6,7]. Thus, coastal wetlands are high-valued ecological zones where terrestrial and marine ecosystems meet and intersect, sharing distinctive features and threats [5,8,9,10]. These ecosystems are often essential refuges for a wide range of flora and migratory coastal birds, as is the case for species that move through the migratory routes of the Pacific coast of South America (Southeastern Pacific Flyway, SEPF) [11,12,13].With nearly 6,500 linear kilometers of continental coastline, the coastal wetlands along continental Chile are thus a crucial element of the SEPF - even with a fragmented distribution that decreases from south to north – since it provides a valuable migratory corridor with high biodiversity, productivity, and endemism [3,4,14].

Despite their importance, coastal wetlands are frequently degraded by a variety of anthropogenic activities including wetland reclamation and land use change, water extraction, and nutrient and pollution overload of incoming waters, often causing cascading effects on abiotic variables and processes, as well as the structure and function of biological hierarchical levels such as primary producers and their consumers (e.g. [15,16,17]. In coastal wetlands, these changes frequently occur in tandem with the effects of global-scale environmental stressors such as the effects of climate change through increasing temperature and sea level rise [5] or the intensity and frequency of storm surges [18]. In Chile, the most vulnerable coastal wetlands are found between 32°S and 40°S, and the dominant drivers affecting them are climate change as well as urbanization and land use changes which affect 41.9% and 52.8% of these wetlands respectively [7]. Among those, the Rio Cruces Wetland (RCW) near the city of Valdivia (ca. 39.9°S stands out, as it harbors Chile’s first Ramsar Site, located in the wetland’s central area. The RCW was formed as a result of large co-seismic continental subsidence following the great Valdivia earthquake of 1960, which caused level decreases of up to 2.5 m [19,20,21]. As a result, land areas that were previously adjacent to the riverine system and occupied by agriculture, stock farming, and marsh forests were transformed into shallow-water or frequently flooded wetland areas [22] (Manzano-Castillo et al. 2020). In 1981, the central area of the RCW (Figure 1) was declared as Ramsar Site n°222, given its support of numerous waterbirds, including two endangered species Coscoroba swan (Coscoroba coscoroba) and the White-faced Ibis (Plegadis chihi). At that time, it harbored the largest recorded stable nesting population of the Black-necked Swan (Cygnus melancoryphus), totaling 3,000 individuals. The RCW is modulated by the tidal cycle, as the water level is affected by inflow of tides through the river mouth at Bahia Corral, with the tidal cycle modulating the water level across the wetland [23]. The RCW currently hosts a high diversity of birds, amphibians, fish, mammals, and aquatic macrophytes [24,25,26].

The RCW has been the subject of several studies following the changes in water quality, after a pulp mill plant (located 15 km upstream the wetland), started operating and discharging residual waters into the river during February 2004 [27]. The changes in water quality cascaded throughout the widespread stands of the Brazilian elodea, Egeria densa (Planchon, 1849; Hydrocharitaceae). This invasive alien species was and remains the main food of the black necked swans at the RCW, and its collapse during 2004 led to massive swan migration and mortality due to emaciation and histopathological liver damage [17,27,28,29,30]. Several studies conducted following the 2004 population decline and subsequent recovery of C. melancoryphus have demonstrated the significance of aquatic macrophytes for the RCW ecosystem and its tributary rivers, as many of them can be considered ecosystem engineers [31,32]. Among these, Egeria densa still stands out as a dominant food source for the herbivorous aquatic birds of the wetland, such as black-necked swans and coots, while it also helps to stabilize the benthic sediments in the wetland’s shallow channels and flatlands, playing important roles in the flow of energy and matter, as well as the structure and hydrological dynamics of this ecosystem [28,31,33]. Thus, this invasive alien species (IAS) plays crucial roles in the RCW ecosystem. Along with this IAS, the native macrophytes Potamogeton lucens (Potamogetonaceae) and Potamogeton pusillus L. (Potamogetonaceae) have also been shown to be important trophic items for swans [33] and stabilizers of sub-aquatic sediments [24,34]. In addition to these important food sources for herbivorous birds and ecosystem engineers, most sedimentary tidal flats in the RCW are associated with or surrounded by the cattail Schoenoplectus californicus (C.A. Mey.) Soják, Cyperaceae) which is one of the RCW’s dominant macrophytes and occupies extensive areas along their subaerial ebb tide perimeter [22]. Today, most of these shallow cattail tidal flat complexes are mostly found in the southern and middle sectors of the RCW [22]. Schoenoplectus californicus, like Egeria densa, stabilizes and models wetland shores or shallow bottoms against external hydraulic stress factors, since its emerged stems or canopy structures capture suspended sediments (cf. [35,36,37,38]. In this regard and given the relevance of the effects of climate change on Chile’s coastal wetlands, it is important to determine the potential effects of climatic and hydrological drivers on these aquatic ecosystems.

Egeria densa is a submerged perennial monocotyledonous dioecious plant indigenous Brazil, Uruguay, and Argentina in eastern South America, which is currently recognized as a problematic invasive alien species by the International Union for the Conservation of Nature [39]. Its popularity as an aquarium species has led to the establishment of wild populations in all the continents, particularly across many temperate regions around the world, including Chile, North America, Oceania, Europe, Africa and Asia [31,39,40,41,42,43,44,45,46,47,48]. As mentioned above, it became the main food source for black necked swans and other herbivorous bird species in the RCW [28,29,30]. Thus, it has become an integral part of the food web of this Ramsar Site and adjacent wetland areas, and in so doing stabilizes shallow bottom sediments in the RCW. On the other hand, Schoenoplectus californicus is a perennial or rhizomatous geophyte. It is a native species widely distributed throughout the Americas and parts of the South Pacific, growing and primarily in the subtropical biome [49]. While it is not an important trophic item for herbivorous waterbirds, it is an important ecosystem engineer and provides habitats both for aquatic species as well as for reed birds, herons and Nutria or Coypu (Myocastor coypus). The contrasting life histories and morphologies shown by both species raise the question as to how they may respond to climatic and hydrologic alterations caused by climate change, and the implication of these responses for the ecosystem engineering functions provided in the RCW by these dominant aquatic macrophytes.

In this regard, the RCW is subject to disturbance effects resulting from anthropogenic factors, including land use change and climate change, as shown previously by [7]. However, other agents of change include alterations in the Earth’s crust due to seismic cycle dynamics [50] and variability impacts from the warm and cold phases of the El Niño-Southern Oscillation (ENSO) phenomenon, that can lead to water deficits in the contributing basins that supply water flow to coastal wetlands [51,52]. Additional sources of climate variation include the variability linked to the Pacific Decadal Oscillation (PDO), which may result in decadal water deficits, as well as the Atlantic Multidecadal Oscillation (AMO), which similarly influences low-frequency variability [53,54,55]. In addition, the Southern Annular Mode (SAM) or Antarctic Oscillation (AAO), also contributes to large-scale climate phenomena by characterizing the north-south movement of the westerly wind belt encircling Antarctica, a phenomenon which significantly influences climate variability in the mid- and high latitudes of the Southern Hemisphere [51,56,57,58]. Those synergic effects between climate change and large-scale climatic forcing have affected different aspects of the RCW’s climate [51,58] (Figure 2).

Considering the differing morphological strategies and life forms of the two dominant aquatic macrophytes of the RCW, the exotic Egeria densa, and the native Cattail, Schoenoplectus californicus, along with their roles as ecosystem engineers in this important Ramsar site and Nature Sanctuary, it is essential to evaluate their responses to climatic and hydrological variations in the wetland area. This study aims to evaluate the relative impact of climatic and hydrological factors on the spatial and temporal distribution patterns of Egeria densa and Schoenoplectus californicus over the last decade. We integrate long term monitoring efforts across a ten-year period (2015 – 2024), including systematic occurrence of data collected during the summer seasons and integrated with remote sensing images, to fit a species distribution model (SDM) for both species in every year. These SDMs allow us to estimate the distribution area occupied by Egeria densa and Schoenoplectus californicus within the Ramsar site, and to model its response to climatic and hydrological variables over the study area. This will allow us to determine what are the relative effects of climatic (temperature and precipitation) and hydrological (flow) drivers on the spatial and temporal distribution patterns of the dominant aquatic macrophytes in the Rio Cruces wetland. By examining these species, this research aims to provide insight into how native and exotic aquatic macrophytes respond to climatic and hydrologic forcing, contributing to a broader understanding of ecosystem dynamics in response to environmental change.

2. Results

2.1. Environmental Variability

The study of the pattern of climatic and hydrologic variability at the study area, showed a clear seasonal pattern in all four environmental variables, with autumn and winter months showing decreases in temperature and increases in precipitation, river flow and water level (Figure 2).

Examination of the general additive models shows that the baseline decades differ significantly from 0 in all four variables (Table 1). The fitted models show significant decreases over the last 10 years (2013 to 2023) in temperature as well as river flow and water level, while precipitation shows a non-significant trend towards earlier rainy seasons (Figure 2, Table 1). The fitted general additive models for temperature, river flow and water level account for a 70-89% of the observed variability, while the seasonal model for precipitation accounts for 55% of the observed variability.

While the observed long-term patterns of change in precipitation for the study area reported previously [51,58] are not reflected in recent changes of the seasonal pattern, the rest of the variables examined do present significant decreases over the recent decade (Figure 2, Table 1). Thus, for these three variables, the models effectively capture the seasonal temperature pattern and the differences between groups. The significant smooth term for month confirms the presence of a strong seasonal trend, while the significant group differences highlight systematic variations between the recent decade versus the baseline. The fact that three out the four climatic and hydrologic drivers examined have changed relative to their long-term baselines suggest these drivers and related variables may be affecting the RCW and its surrounding hydrologic basin, with potential effects on the dominant aquatic macrophyte species of the wetland.

2.2. Species Distribution Modelling

Despite the high frequency of cloud cover across the study area, it was possible to obtain cloud-free Landsat 8 scenes in all ten spring-summer study seasons (Table 2).

Across all ten years and for both species, the fitted SDM models showed very good fits to the observed presence data, with Area Under the Receiver Operating Characteristic Curve (AUC) values above 0.80 both the training and test cross-validation data sets across all samples for Egeria densa and Schoenoplectus californicus presenting AUC values above 0.85 (Table 1), which indicates the models provide very good fit in the classification task, while those years with test AUC values greater than 0.9 present highly reliable fits to the available data [59]. Detailed Receiver Operating Characteristic Curve (ROC) curves for Egeria densa and Schoenoplectus californicus are shown in Supplementary Figures S1 and S2 respectively. The spatial pattern of the fitted MaxEnt habitat suitability maps is shown in Figure 3 and Figure 4. In the case of the exotic macrophyte Egeria densa, we observed a heterogeneous spatial variation in habitat suitability (Figure 3). Thus, patches of high habitat suitability across all the shallow areas of the wetland were observed during the first four years (2015 to 2018). After 2019, habitat suitability decreases in the northern sections of Rio Cruces and northern tributaries, followed by a marked decrease in habitat suitability across the central and southern areas of the wetland (Figure 3e–g). Some of the areas that show decreases in habitat suitability do show signs of increase after 2021, suggesting a fluctuation in the distribution of this species (Figure 3h–j). On the other hand, results for the native the cattail, Schoenoplectus californicus show a patchy distribution across the wetland, associated with shallow tidal flats, with decreases in the values of habitat suitability in the north section of the wetland being observed in the spring-summer seasons of 2017, 2020, 2021 and 2022 (Figure 4a–h). The fluctuation in habitat suitability in northern sections of the wetland contrast with the areas in the central and southern parts of the wetland, which tend to present stable patches of the cattail (Figure 4).

2.3. Assessing Interdecadal Variation in the Area of Suitable Habitat

We then examined the dynamics of the area of suitable habitat within the RCW. For both Egeria densa and Schoenoplectus californicus, we observed that while the area shows yearly oscillations, no significant linear trend is present (Figure 5). Thus, the Ordinary least Squares (OLS) linear regression for Egeria densa was not statistically significant (R² = 0, F(1,8) = 0.005, p = 0.943), with no significant effect of year on area (β_Year = -4±53.8, p = 0.943). In the case of Schoenoplectus californicus, while a decreasing trend was observed, this was not statistically significant (R² = 0.22, F(1,8) = 2.3, p = 0.168), and again, no significant effect of year on area was observed (β_Year = -87.9±58, p = 0.168) (Figure 5). When examining the degree of relative variability Egeria densa showed a lower Coefficient of Variation (CV) value than Schoenoplectus californicus (20.46% and 30.05% respectively), suggesting that these species show different responses to environmental fluctuations.

After the above-mentioned analyses, we conducted a stepwise linear regression analysis to evaluate the interannual variation in macrophyte area within the RCW as a function of climatic and hydrological variables. For the IAS Egeria densa, the final model retained three significant predictors: mean annual temperature (T), standard deviation of water level (s.d. Level), and Year. The model intercept was significant, -3733.32 ± 1157.59 (t = -3.23, p = 0.018) (Table 3). In this model, the Year had a positive association with macrophyte area (β = 1.83 ± 0.57, t = 3.22, p = 0.018), indicating an increasing trend over time. Mean annual temperature (T_Year) also showed a significant positive effect (β = 6.93 ± 1.72, t = 4.02, p = 0.007), suggesting that higher temperatures contribute to an expansion of Egeria densa. In contrast, standard deviation of water level (s.d. Level) had a significant negative effect (β = -32.94 ± 11.76, t = -2.80, p = 0.031), indicating that greater fluctuations in water level may negatively impact the species’ area of suitable habitat.

For the native Schoenoplectus californicus, the stepwise linear regression model retained all five predictors: year, mean annual temperature (T_Year), accumulated annual precipitation (sP_Year), mean annual water level (Level), and standard deviation of water level (s.d. Level). While the intercept and year were not statistically significant (p = 0.093 and p = 0.098 respectively), mean annual temperature (T_Year) had a significant negative effect on Schoenoplectus californicus (β = -6.19 ± 1.85, t = -3.35, p = 0.029), suggesting that higher temperatures reduce the area of suitable habitat. Accumulated annual precipitation (sP_Year) was positively associated with macrophyte area (β = 0.02 ± 0.003, t = 5.13, p = 0.007), indicating that increased rainfall favors the expansion of Schoenoplectus californicus. However, mean annual water level (Level_Year) had a significant negative effect (β = -46.47 ± 9.99, t = -4.65, p = 0.010), implying that higher river levels are associated with a reduction of suitable habitat area for Schoenoplectus californicus. Lastly, standard deviation of water level (s.d. Level) was not a significant predictor (p = 0.294) (Table 3).

For Schoenoplectus californicus, the selected regression model included non-significant effects (see Table 3). Hence, we also examined a simplified linear model was fitted that excluded those non-significant variables. This enabled us to compare the stepwise linear regression model to the reduced model, with the Akaike Information Criterion (AIC) determining and selecting the most parsimonious model [60,61]. The reduced model exhibited a higher Akaike Information Criterion (AIC = 55.4) compared to the stepwise model (AIC = 51.7), suggesting that the inclusion of all five predictors provided a better overall model fit, while remaining parsimonious (see supplementary Table S1 in Supplementary Materials). Despite this, both models identified similar key drivers influencing Schoenoplectus californicus habitat. Mean annual temperature (T_Year) consistently had a significant negative effect, indicating that higher temperatures reduce suitable habitat area. Similarly, accumulated annual precipitation (sP_Year) had a strong positive influence, while mean annual water level (Level_Year) had a significant negative effect, suggesting that increased water levels reduce the area of suitable habitat. These results confirm that climatic and hydrological variables play a crucial role in determining the extent of Schoenoplectus californicus habitat, and that while the reduced model simplifies interpretation, the stepwise model provides a more comprehensive representation of environmental influences.

Overall, the analysis reveals that the IAS Egeria densa responds positively to increasing temperatures and tends to increase over time, while native Schoenoplectus californicus is negatively affected by higher temperatures and river levels but benefits from increased precipitation. The significant predictors identified in these stepwise regression models highlight the differential responses of these two macrophyte species to climatic and hydrological variations.

3. Discussion

3.1. Climatic and Hydrological Variability in the Study Area

Our results show significant seasonal and interannual climatic and hydrological variability within the RCW, with distinct seasonal patterns in air temperature, precipitation, river flow, and water level. While average monthly air temperature, river flow, and water level have all decreased significantly over the last decade (2013-2023), precipitation does not show significant changes, although a trend towards earlier rainy seasons can be appreciated. These findings are consistent with earlier research showing long-term alterations in precipitation patterns and hydrological regimes in the region [51,52,56,57,58]. On the other hand, the observed reduction in air temperature contradicts global and regional warming trends [62], which suggests the possible effects of localized climate influences on the study area, likely resulting from the influence of neighboring oceanic water mass, or interactions with meteorological events such as atmospheric rivers [63,64,65]. While recent studies have examined the effects of climate change and hydrological variation on watersheds across Chile (e.g., [66]), the detailed variation of climate and hydrology in some coastal watersheds, such as the section studied of the RCW, have yet to be investigated. Our findings highlight the need for additional research into the impact of climate change on climatic and hydrologic variables in this area of the RCW. This is particularly relevant in the light of the recent accreditation of Valdivia as one of Latin America’s first Ramsar wetland cities.

3.2. Species Distribution Modelling

Despite the high cloud cover frequency of the region, cloud-free Landsat 8 scenes were acquired for every study year. The increasing availability of Landsat and Sentinel imagery would increase the availability of these spectral satellite imagery, increasing the likelihood of expanding these modelling efforts to autumn and early spring seasons of the RCW. Thus, our results show how remote sensing and species distribution modelling (SDM) may allow successful monitoring of aquatic macrophytes in similar wetland habitats to those of the RCW [67,68,69].

Across all ten years analyzed, the fitted species distribution models for both the exotic Egeria densa and the native Schoenoplectus californicus demonstrated high predictive performance, with test AUC values exceeding 0.80. Notably, 7 out of 10 Egeria densa MaxEnt models and 8 out of 10 Schoenoplectus californicus models, achieved test AUC values of 0.90 or higher, confirming their robustness in assessing species-environment relationships and distinguishing between suitable and unsuitable habitats (cf. [70,71,72]. The strong predictive capacity of these remote sensing-based SDMs highlights their utility for understanding macrophyte distribution patterns in the RCW, and their potential application in wetland monitoring and conservation planning in other similar ecosystems. While some authors have pointed out some limitations on the use of MaxEnt as a SDM algorithm [73], different large-scale assessments have validated its usefulness and performance to model presence-only data with small data sets [71,72,74].

3.3. Assessing Interdecadal Variation in the Area of Suitable Habitat

Egeria densa experienced a decline in habitat suitability after 2019, followed by some recovery post-2021, while Schoenoplectus californicus showed stability in the central and southern wetland but fluctuations in the northern sections. As a result, no significant linear trend in was detected for either species, although both present marked interannual variability, with Egeria densa showing lower relative variability than Schoenoplectus californicus. However, stepwise regression analyses revealed that these two species differ in their responses to climatic and hydrological factors. For Egeria densa, higher mean annual temperature emerges as a likely driver of suitable habitat expansion, while greater water level fluctuations negatively impacted its distribution. These effects are consistent with previously documented stress response patterns Egeria densa in Japan [45,46,47]. By measuring the production of hydrogen peroxide (H₂O₂) in plant tissues as a marker for stress response, with plants of Egeria densa, showing reduced growth and a deteriorated physiological condition when tissue H₂O₂ concentration exceeds a threshold value [45,46,47]. Field measures across different rivers showed that for temperatures between 10 to 30°C, H₂O₂ concentration in Egeria densa tissues decreased with increasing temperature [46]. Furthermore, H₂O₂ concentration increased with turbulent flow velocity [46]. Increases in water level fluctuations may be associated both with increased flow variability and average value (See Supplementary Figure S4), as well as with more extreme turbulent flow events. Both conditions could account for the observed negative effect on suitable habit area at the wetland level. Thus, our observed results for Egeria densa are consistent with available information on the stress drivers for this exotic species. These results and previous observations suggest that other wetlands and continental water areas with low turbulent flow, such as lakes or rivers and channels with low slopes may have greater habitat suitability for this species. Furthermore, any future changes in water temperature patterns may also alter the habitat suitability for this exotic species. Further studies are needed to determine whether the pattern of H₂O₂ concentration previously described in [45,46,47] can also be observed in the RCW wetland or other invaded water courses in South America.

In contrast with the invasive Egeria densa, native Schoenoplectus californicus, exhibited a negative association with temperature and mean annual water level, along with a positive association with accumulated precipitation. This suggests that decreasing temperatures and mean annual water level, may increase suitable habitat for this native species. The observed variation - in both temperature and mean annual water level - suggests that interannual variability in both environmental drivers may account for observed changes in the available area of suitable habitat. Our results are partially consistent with experimental studies conducted in the Sacramento - San Joaquin Bay Delta, California, which show that while transplanted Schoenoplectus californicus plants can tolerate more severe frequency, depth, and duration of flooding than other aquatic macrophyte species (Schoenoplectus acutus and Typha latifolia), this species showed greater vegetation expansion in transplant sites characterized by a deeper surface layer of non-compacted soil in conjunction with shorter durations of flooding [75]. Thus, decreased mean annual water level may result in shorter durations of flooding in fringe habitats, allowing greater expansion to occur [75]. Furthermore, experimental greenhouse studies of different life history stages of Schoenoplectus californicus and Schoenoplectus acutus, shows that longer flooding durations led to lower survival rates for seedlings of both species. These results would be consistent with the conditions of increased mean annual water level, which would entail longer flooding events and a decrease in seedling survival [76]. However, further studies are required that examine the variation in the degree of soil compaction and flooding dynamics (cf. [76]), in order to know whether these variables are the proximal cause of observed spatial and temporal variation in the cover and survival of Schoenoplectus californicus across the RCW. Our results show that successful management and conservation of Chile’s first Ramsar site will benefit from the integration of long term monitoring efforts with remote sensing and modelling strategies to develop useful insights on the potential consequences of ongoing climatic and hydrological changes at the watershed level.

4. Materials and Methods

4.1. Study Area and Environmental Variability

The Rio Cruces Wetland (RCW) is situated in the coastal area of the Valdivia River basin, and is characterized by a temperate rainy climate, experiencing significant rainfall during the winter months, with limited dry periods and low temperatures in winter [51]. The hydrological dynamics of RCW sub-basin and its contributing river basins are primarily influenced by precipitation [52,66,77]. The RCW and its surrounding hydrographic basin are situated in a temperate macro bioclimate, significantly influenced by the hyper-oceanic temperate bioclimate and characterized by low thermal oscillation [78]. Remarkably, over 50% of the annual precipitation occurs from May to August [51].

4.1.1. Field Surveys

Starting in 2014, field monitoring was conducted in the Rio Cruces Wetland (RCW) to document the presence of dominant aquatic macrophytes, with particular focus on the exotic invasive species, the Brazilian Elodea, Egeria densa Planch. (Luchecillo in spanish, [49] as well on the native Cattail (also known as California bulrush), Schoenoplectus californicus (C.A.Mey.) Soják (Totora in spanish, [49]). These are two contrasting species which play important roles in structuring the RCW ecosystem, as they are either an important food source for herbivorous waterbirds (Egeria densa) or play an important role as ecosystem engineers in the RCW (Egeria densa and Schoenoplectus californicus). Geo-referenced occurrences of these species have been sampled during the austral spring-summer seasons between 2014 and 2023, yielding a total of 10 years. For each of these years, we located and recorded large patches of these two species in the RCW, selecting patches with diameter equal or greater than 30 m, which corresponds to the spatial resolution of Landsat scenes. For each of these stands, geographic coordinates of the center of the stand were registered using a Global Positioning System navigator (GPS), with a WSG84 coordinate system with a UTM datum (18S zone). The presence data for large monodominant patches were then used to fit species distribution models for each species across the wetland in each of these 10 spring-summer seasons.

4.1.2. Environmental Variability

To describe the climatic variability in the Rio Cruces Wetland (RCW), we obtained and analyzed hourly time series records of air temperature (T_hour, measured in °C) and precipitation (P_hour, measured in mm) from the Pichoy Airport meteorological station (39.65667°S, 73.08722°W), located 22 km northeast of the city of Valdivia (Figure 1). To analyze the hydrologic variability in the area where the RCW is located, we analyzed daily time series of river flow (m³/s) and water level (m) data from the Rucaco hydrological station (39.55°S, 72.90°W), located 42 km northeast of Valdivia, in the northern upstream reaches of the Rio Cruces River (Figure 1). Air temperature and precipitation records at Pichoy Airport are available from January 1st, 1966, while daily river flow and water level at Rucaco hydrologic station are available from May 1st, 1969, and January 1st, 2000, respectively. Time series were analyzed using records up to December 31st, 2023. Each hourly time series (T and P) was inspected to identify regions with missing values, which were then interpolated by fitting a structural time series model to capture observable dynamics and impute missing data. Missing value imputation was carried out using R’s imputeTS library’s na_kalman function was used for this approach [79,80]. This involves applying a Kalman filter to represent the time series and interpolating missing observations. The Kalman filter is a technique for estimating the dynamics of a linear dynamical system using partial observations with white additive noise or equal variance at all frequencies. The former is useful when noise or inaccuracy is detected in measurements and an accurate estimate of the system’s dynamics is required [81]. Once the missing values were imputed in the hourly T and P time series, these were summarized at a daily level, calculating mean daily values. These daily time series for T, P, Flow and Level were then used to (i) assess recent interdecadal variation relative to the long-term available baseline, and (ii) assess recent patterns of variability of climatic and hydrological drivers across the studied ten-year period.

4.1.3. Assessing Interdecadal Variation

To assess the interdecadal variation, daily time series of all four variables (T, P, Flow and Level) were summarized to monthly average values. These monthly time series were then modelled by using Generalized Additive Model (GAM) framework with cyclic cubic splines, which are useful for periodic data (e.g., months in a year, hours in a day). This ensures that the start (January) and end (December) of the cycle match smoothly, describing any phenological or seasonal cycles present. We included a Group linear predictor which identified data corresponding to the baseline and data in the most recent decade (2013–2023). The baseline was defined by considering the starting date for each time series. All analytic procedures were carried out using R, with mgcv and ggplot libraries [80,82,83,84,85,86].

Assessing Recent Climatic and Hydrological Variability

To examine recent climatic and hydrological variability in climatic and hydrological drivers that may affect wetland structure and dynamics, we examined the monthly time series and calculated a set of potentially relevant variables. Given that the field sampling was done on late spring or summer austral seasons (see Table 1), we focused on data for the years prior to the studied sampling dates. Thus, climatic and hydrologic variables were calculated for the years 2014 to 2023. Based on the monthly time series, we calculated mean annual values for all four variables (T, P, Flow and Level), as well as annual standard deviation values for P, Flow and Level. In addition, we also calculated accumulated annual precipitation values as well as the year as a time variable. The resulting nine annual time series were then examined to exclude any variables with high degree of multicollinearity, measured as the variance inflation factor [87,88], by using a threshold value of 10, as commonly used in the literature. This procedure was implemented using the vifstep function in R library usdm [88]. As a result, five variables were retained: Year (VIF = 6.05), mean annual Temperature (VIF = 5.77), mean annual water level (VIF = 3.64), annual standard deviation of water level (VIF = 2.30) and cumulative annual precipitation (VIF = 5.18). All analytic procedures were carried out using R [80].

4.2. Remote Sensing Image Acquisition and Processing

To fit SDMs for each species in each of the ten austral spring-summer sampling seasons, we used remote sensing layers extracted from a Landsat 8 operational land imager (OLI) scene. In each spring-summer sampling season, we downloaded the closest available OLI scene located in path/row combination 233/88 of the Worldwide Reference System 2 (WRS-2). This scene encompasses the entire RCW and is centered at 40°19’20” S, 72°51’00” W. We processed bands 2–7 of the OLI scene to generate a set of SDM predictor Geographic Information System (GIS) layers, following previous studies in this wetland [17,28,29]. Briefly, bands were radiometrically calibrated using Landsat 8 radiance rescaling factors provided in the Landsat image metadata file, and top of atmosphere spectral radiance values (L_λ, W·(m2·sr·μ m)-1) were calculated. These L_λ values were then converted to top-of-atmosphere reflectance percentages (R_TOA) [89], and the atmospheric correction for Case-2 turbid waters was applied using the path extraction method [17,90]. Following these corrections, the bands were clipped to the study area, and two spectral indices were calculated for each scene: a chlorophyll proxy (CHL = Blue/Green = Band 2/Band 3) [17], and the normalized difference vegetation index (NDVI = (NIR - Red)/ (NIR + Red) = (Band5 − Band4) / (Band5 + Band4), where NIR corresponds to the Near Infrared band [91]. To limit the modeling to wetland areas and river courses, a river raster layer that masks terrestrial pixels was generated using existing cartography of water courses in the Region, as well as previous raster layers generated by previous Landsat-based studies [92,93,94].This resulted in a set of nine predictive GIS layers with a spatial resolution of 30 m pixel size that characterize the studied scenes. For each of the ten years studied, all nine GIS layers were used to fit a SDM that provides a spatially explicit estimate of macrophyte habitat suitability distribution, as described in the following section. All GIS layer processing was carried out using. GIS and statistical data processing were conducted utilizing the dplyr, terra, and ggplot2 packages within the R statistical computing environment [86,95,96] as well as QGIS version 3.20.3 Odense [97].

4.3. Species Distribution Modelling

SDMs were fitted using Maximum Entropy species distribution modelling software (MaxEnt v.3.3). MaxEnt uses information on spatial occurrences or presences and GIS layers or features to estimate the presence probability function across the study area [71,72,74,98,99,100,101]. Comparative studies have shown that MaxEnt performs better in relation to other methods developed to analyze presence only data sets (e.g., [72,74,102], even in situations where the sample size is small [103,104,105]. The resulting maximum entropy model has been shown to be equivalent to maximizing the likelihood function of a spatial inhomogeneous Poisson point process [106,107,108], and as a result, its output can be interpreted as providing relative density estimates across space [106,108]. Hence, the resulting probability measure provides a continuous estimate of probability, which measures environmental habitat suitability (HS) for the species under study. MaxEnt SDMs for each of these two macrophytes in the 10 years studied were fitted using a 5-fold cross-validation scheme, thus allowing every occurrence data point to be used as part of the training and evaluation data sets [71,72,74,98,99,101]. To measure how well each model discriminates presences more accurately than a random prediction, we used the area under the curve (AUC) statistic for the Receiver Operating Characteristic (ROC), with values greater than or equal to 0.8 being considered indicative of a very good or excellent fit, while values greater than 0.9 are considered as evidence of an outstanding or highly reliable fit to the available data [59]. Fitted models were later projected over the Rio Cruces wetland, using the same GIS predictive layers, with the resulting species distribution HS map allowing us to identify where suitable environments for each of the aquatic macrophytes are found in the wetland.

4.4. Statistical Analysis

The overall variation in distribution among the two study species was summarized by analyzing the spatially explicit estimation of probability (or HSI values) across the wetland. First, we generated visual representations of the time series of probability maps for each aquatic macrophyte species. Second, we estimated the suitable area available in each year for both the exotic Egeria densa and the native Schoenoplectus californicus. To do this, each annual probability map for each species was converted into a presence-absence map by applying a threshold value τ to the predicted presence probabilities of each species. This threshold was chosen by using a probability threshold value that maximizes the sum of sensitivity and specificity (MSS) [109]. The resulting time series of presence-absence maps delineates the estimated distribution dynamics for both macrophytes. We then estimated the area of available suitable habitat for each species over the study period, by calculating the number of occupied pixels in each year and subsequently multiplying this figure by the surface area of a pixel in a Landsat image (900 m²). This yielded a time series of estimated areas of suitable habitat for each species.

To quantify the degree of interannual variability in area of suitable macrophyte habitat over time within the RCW, we calculated the coefficient of variation (CV), defined the ratio of the standard deviation to the mean, expressed as a percentage. This allowed us to assess the relative fluctuations in habitat availability across the two macrophytes studied [110]. A low CV indicated stability in suitable habitat area, suggesting consistent environmental conditions, whereas a high CV reflected significant temporal variation, potentially driven by climate variability, hydrological cycle changes, or other ecological disturbances. This approach allowed for a standardized comparison of habitat variability across different regions and time periods. To assess the relative effect of climatic and hydrological variables, we then conducted a stepwise linear regression analysis separately for Egeria densa and Schoenoplectus californicus. In each species’ regression model, the response variable was the annual area of suitable habitat, while the predictor variables were selected after removing any collinear variables (see section 4.1.2.2). The selected variables were Year, mean annual temperature (T_Year, °C), accumulated annual precipitation (sP_Year, mm), mean annual water level (Level_Year, m), and annual water level standard deviation (s.d. Level, m). Stepwise regression was implemented using the function stepAIC from R’s MASS library [111]. This function allows us to iteratively select the most relevant predictors based on statistical significance, optimizing model fit while minimizing multicollinearity. We then examine the plots of observed and predicted annual area of suitable habitat, comparing it with the linear model prediction based on Year alone. In those cases where the selected regression model included non-significant effects, a reduced linear model was fit, excluding non-significant variables. This allowed us to compare the stepwise linear regression model with the reduced model, using the Akaike Information Criterion (AIC) to identify and select the most parsimonious model [60,61].