Rapid plant functional trait responses to warming, flooding, and herbivory in high-latitude coastal wetlands

TL;DR

This study shows how plants in high-latitude wetlands quickly adjust their traits in response to climate change factors like warming, flooding, and herbivory.

Contribution

The study reveals context-dependent and consistent trait responses to multiple climate drivers in high-latitude wetlands using a full-factorial experiment.

Findings

Warming increased plant size in the Lowland wetland but not in the Upland wetland.

Herbivory consistently reduced plant size and promoted acquisitive leaf traits across all species.

Flooding had contrasting effects on leaf economics traits depending on wetland type and species.

Abstract

Climate change is rapidly altering high-latitude coastal wetlands through increasing temperatures, more frequent flooding, and changing herbivore abundance and distribution. Examining plant functional trait responses to these drivers provides insight into how plant communities are shaped by climate change. We used a one-year full-factorial mesocosm experiment in coastal Alaska to assess immediate responses of two size-related traits (vegetative height and leaf area) and two leaf economics traits (specific leaf area [SLA] and leaf dry matter content [LDMC]) to warming, flooding, and goose herbivory in three dominant species across two wetland communities (Lowland and Upland) differing in landscape position. We sampled different sedges (Carex rariflora and Carex lyngbyei) and the same deciduous dwarf-shrub (Salix fuscescens) in the two wetlands. Warming increased plant size for the sedge (leaf area: + 15%) and the shrub (height: + 15%; leaf area: + 19%) in the Lowland wetland only. Flooding promoted acquisitive trait values for the Lowland sedge (SLA: + 8%; LDMC: − 6%) and conservative values for the Upland shrub (SLA: − 9%; LDMC: + 11%). The shrub thus showed clear wetland-specific responses to warming and flooding. Herbivory triggered coordinated responses across all three species, consistently decreasing size-related trait values (height: − 17–35%; leaf area: − 33–50%) and promoting acquisitive trait values (SLA: + 14–40%; LDMC: − 12–23%). Despite the absence of interactions, each driver contributed significantly to trait variation, underscoring the value of multifactorial approaches. The combination of consistent (herbivory) and context-dependent (flooding and warming) responses highlights the complexity of trait responses and improves predictions of rapid phenotypic adjustments in coastal high-latitude wetlands.

The online version contains supplementary material available at 10.1007/s00442-026-05876-8.

Figures

Plant species in the multivariate functional trait space. PCA biplot of plant traits: height, leaf area (LA), specific leaf area (SLA), and leaf dry matter content (LDMC), shown as black arrows. Each point represents the PCA scores of one species in a mesocosm and is colored by species identity. *S. fuscescens* is shown separately for the two wetlands. Ellipses indicate 95% confidence intervals around species centroids. Percent variance explained by the first two principal components is shown on the axes

Effects of treatments on multivariate functional trait expression. PCA biplots and PERMANOVA results of plant traits under warming, flooding (low- and high-intensity), and herbivory. Panels show results for: (A) *C. rariflora*, (B) *S. fuscescens* in the Lowland wetland, (C) *C. lyngbyei*, and (D) *S. fuscescens* in the Upland wetland. While treatment-specific biplots are displayed separately for interpretability, PCA and PERMANOVA were conducted once per species. Ellipses indicate 95% confidence intervals around treatment centroids. The four focal traits are shown as black arrows in the left

Two-way interaction between flooding and warming on vegetative height of *S. fuscescens* in the Upland wetland. GLMM predictions (± 95% confidence intervals) under combinations of flooding and temperature treatments. Treatment levels include ambient temperature (A), warming (W), no flooding (nF), low-intensity flooding (LF), and high-intensity flooding (HF). Bolded values indicate marginally significant effects (*P* < 0.10); asterisks denote significance at *P* < 0.05

Full text

Introduction

Coastal high-latitude regions are undergoing rapid climate change (Irrgang et al. 2022; IPCC 2023). The Arctic as a whole has warmed nearly four times faster than the global average (Rantanen et al. 2022). Concurrently, coastal flooding has intensified (Kirezci et al. 2020) due to accelerated sea-level rise (Hamlington et al. 2024), diminishing sea ice extent and increased storm activity (Vermaire et al. 2013), land subsidence (Nicholls et al. 2021), and more frequent, intense cyclones (Parker et al. 2022). At the same time, goose herbivore pressure is shifting as warming advances migration phenology, overall abundances are increasing, and rising coastal flooding is likely driving geese toward previously less utilized habitats further inland (Fox and Madsen 2017; Lameris et al. 2021; Koltz et al. 2022). Understanding how these concurrent global change drivers influence plant communities in high-latitude coastal ecosystems is critical for anticipating ecological change in these rapidly transforming regions.

A deeper understanding of climate change impacts on plant communities can be gained by examining how functional traits respond to environmental change (Suding et al. 2008; Myers‐Smith et al. 2019), as these traits mediate species’ interactions with their environment and determine their capacity to adjust to new conditions (Violle et al. 2007). Plants with greater trait plasticity are generally better able to cope with changing environments, while those with limited capacity for phenotypic adjustment may be more vulnerable (Nicotra et al. 2010; Henn et al. 2018). Aboveground functional trait variation is primarily structured along two axes: plant size and the leaf economics spectrum (Díaz et al. 2016). Size-related traits such as vegetative height and leaf area are key to light competition and reproductive success (Pérez-Harguindeguy et al. 2013; Beccari and Carmona 2024). In contrast, the leaf economics axis reflects a trade-off between rapid resource acquisition and efficient resource conservation (Wright et al. 2004; Reich 2014). It is commonly captured by specific leaf area (SLA) and leaf dry matter content (LDMC), with high SLA and low LDMC indicating acquisitive strategies (fast growth, high photosynthetic rates) and the opposite combination indicating conservative strategies (resource retention, stress tolerance). Investigating how these trait dimensions respond to warming, flooding, and herbivory in high-latitude coastal plant communities may help predict functional and ecosystem-level consequences of ongoing environmental change.

Low temperatures make Arctic ecosystems generally temperature-limited, such that even modest warming can alleviate cold constraints and elicit substantial changes in plant growth and associated traits. Warming effects on plant traits in the Arctic are generally well-documented (Myers‐Smith et al. 2019), with consistent evidence for increased height and leaf area (Hudson et al. 2011; Baruah et al. 2017; Bjorkman et al. 2018). Reduced temperature limitation also tends to favor a shift in economics traits toward resource acquisition (Bjorkman et al. 2018; Jessen et al. 2020; Wei et al. 2023), although more conservative responses have also been observed (Hudson et al. 2011). In fact, responses in both size-related and especially economics traits vary with environmental conditions and species identity (Baruah et al. 2017; Betway et al. 2021; Kemppinen and Niittynen 2022). Across the tundra biome, shrubs generally exhibit stronger growth responses to warming than graminoids, particularly in wetter and warmer sites (Elmendorf et al. 2012), and the positive relationship between SLA and temperature is also stronger in wetter environments (Bjorkman et al. 2018). These findings underscore the importance of considering environmental context and co-occurring global change drivers, such as flooding, which might modify the conditions under which warming influences functional trait responses in coastal high-latitude wetlands.

Increased flooding can affect plant functional traits both indirectly, by altering environmental conditions, and directly through physical stress. Depending on whether the benefits (e.g., higher moisture and nutrient availability) outweigh the stresses (e.g., oxygen deprivation, sedimentation, salinity, and mechanical disturbance), flooding may induce contrasting responses in size-related and leaf economic traits. Accordingly, trait responses to flooding have varied across ecosystems, with plant height shown to decrease (Baastrup-Spohr et al. 2015; Fu et al. 2015) or increase (Purcell et al. 2019), and SLA and LDMC either increasing (Jung et al. 2010; Violle et al. 2011) or decreasing (Fu et al. 2015; Purcell et al. 2019). Such divergent outcomes may reflect species-specific differences in flood tolerance. For instance, in high-latitude coastal ecosystems, flooding has been shown to reduce shrub biomass while enhancing graminoid growth (Person and Ruess 2003; Churchill et al. 2015; Petit Bon et al. 2024), suggesting that flood-tolerant, fast-growing graminoids may adopt more acquisitive trait expressions under increased flooding, whereas slower-growing, less flood-tolerant shrubs may display reduced growth and more conservative trait values. As flooding effects on plant traits remain poorly studied in high-latitude coastal ecosystems, further investigation is essential for understanding functional responses to this intensifying global change driver.

In addition to warming and flooding, herbivores can shape plant functional traits by removing biomass through grazing, trampling vegetation and soil, and returning nutrients via feces and urine. Herbivory generally reduces size-related trait values through direct biomass loss, disturbance from trampling, and the selection of larger plants (Díaz et al. 2007; Egelkraut et al. 2020; Barbero-Palacios et al. 2024). Moreover, avoidance-based species—those minimizing the likelihood or intensity of damage—typically display more conservative traits (low SLA and high LDMC). In contrast, grazing-tolerant plants—those capable of compensating or regrowing after tissue loss—tend to express acquisitive traits that support rapid regrowth, often enhanced by nutrient inputs from animal excreta (Díaz et al. 2007; Niu et al. 2016; Jessen et al. 2020). Graminoids, being more grazing-tolerant and capable of rapid regrowth, might thus be expected to shift toward more acquisitive traits, whereas shrubs, which are less able to replace lost tissue, may respond with more conservative trait shifts that reduce further damage. Yet, while geese promote regrowth in preferred forage species such as graminoids (Cargill and Jefferies 1984; Beard et al. 2019), how they affect their functional traits, and how these responses compare to those of the less preferred shrubs, remain unclear. Given the central ecological role of herbivores in high-latitude ecosystems, assessing their influence on plant traits, especially alongside other global change drivers, is a critical research priority (Barrio et al. 2025).

The Yukon-Kuskokwim (Y-K) Delta in western Alaska, one of the largest high-latitude riverine deltas in North America (~ 130,000 km^2^), is undergoing rapid change. The region is warming rapidly (SNAP 2020) and increasingly exposed to frequent flooding, including prolonged overbank flow from monthly high tides and intensified storm events (Jorgenson & Ely 2001; Terenzi et al. 2014). It also serves as a major summer breeding ground for hundreds of thousands of migratory geese (Lyons et al. 2024), which strongly influence vegetation (Uher‐Koch et al. 2019). These factors make the Y-K Delta an ideal system to investigate how warming, flooding, and herbivory, individually and in combination, shape plant functional trait expression, with findings likely relevant to other high-latitude coastal ecosystems facing similar environmental pressures.

In this study, we examined responses of four key plant functional traits related to size and leaf economics to warming, flooding, and herbivory in the wetlands of the Y-K Delta. Using a one-year mesocosm experiment in two distinct wetland plant communities, we applied a full-factorial manipulation of these global change drivers to assess immediate trait responses, reflecting trait plasticity, in three dominant species (two sedges and a deciduous dwarf-shrub). Because these drivers are expected to co-occur in the landscape, understanding whether their combined effects are additive or interactive is essential for predicting species responses. We focused on dominant species, whose traits disproportionately influence ecosystem functioning (Grime 1998) and selected two different functional groups to capture variation in life history strategies and resource acquisition. The presence of the same shrub species in both wetlands also allowed us to assess intraspecific trait plasticity across different biotic and abiotic conditions. We hypothesized that: (1) warming would generally increase plant size and favor acquisitive trait values, particularly in the deciduous shrub; (2) flooding would increase plant size and favor acquisitive trait values in graminoids but may reduce size and favor conservative traits in the shrub, indicative of stress responses; and (3) herbivory would reduce plant size while promoting acquisitive trait values, especially in graminoids.

Methods

Study area and wetland plant communities

We conducted this study in the Yukon Delta National Wildlife Refuge, a vast wetland system (~ 75,000 km^2^) formed by the Yukon and Kuskokwim rivers at their confluence with the Bering Sea (Appendix 1: Fig. S1). Our study sites were located 19 km inland from the Bering Sea and 45 river-km upstream from the mouth of the Kashunuk River. The climate is cold oceanic, with mean summer (June–August) and winter (January–March) temperatures of 12.5 °C and − 12.2 °C, respectively, and an average annual precipitation of 499 mm (30 year average for 1991–2020 at the Bethel weather station, ~ 200 km from the study site; Palecki et al. 2021). The plant growing season spans from late May to mid-August, while snow and ice typically cover the ground from October to early May.

We conducted the experiment in two distinct wetland plant communities, which we refer to as Lowland wetland (61° 26′ 10.93"N, 165° 26′ 20.98"W; 2.39 m elevation) and Upland wetland (61° 26′ 7.44"N, 165° 26′ 36.10"W; 2.56 m elevation). Both wetlands are characterized by fully saturated, oligohaline soils made of a relatively thick organic layer, and both experience occasional summer flooding. However, due to its slightly lower elevation, the Lowland wetland is subject to more frequent flooding, which contributes to differences in plant community composition and soil properties compared to the Upland wetland.

While both wetlands are dominated by Carex species and include dwarf-shrubs such as the relatively abundant Salix fuscescens, the subdominant Betula nana, and the forb Potentilla palustris, C. rariflora dominates the Lowland wetland, whereas C. lyngbyei dominates the Upland wetland. Soil properties also differ between the two wetlands: organic matter content (58% in Lowland vs. 40% in Upland), soluble salts (0.8 vs. 0.5 mmho cm^−1^), sodium (700 vs. 350 ppm), nitrate (0.18 vs. 0.10 ppm), phosphorus (10 vs. 14 ppm), and sulfate (130 vs. 70 ppm). In contrast, soil pH (~ 5.5) and C:N (~ 22) were similar in both wetlands (Ross 2024). These two communities are representative of the widespread wetland communities that compose the complex plant community mosaic of the Y-K Delta (Jorgenson 2000).

Experimental design and treatments

We conducted the experiment within an 8 m × 16 m area in each of the two wetlands with homogeneous plant community composition and minimal topographic variation. At the end of the 2022 growing season, we established mesocosms by removing 60 cm × 41 cm × 18 cm (length × width × height) soil and vegetation monoliths and placing them into transparent plastic containers of the same dimensions. These containers were reinserted into the ground at the same location from which the monoliths had been extracted (Appendix 1: Fig. S2). Each wetland contained 60 mesocosms, for a total of 120 across both sites. Mesocosms were arranged in five rows of 12, with 1 m spacing in all directions. The setup was conducted the year before measurements to allow the monoliths time to adjust to the experimental conditions. Treatment implementation and data collection took place in summer 2023.

Within each row, treated separately to avoid spatial clustering of treatment combinations by chance, we randomly assigned mesocosms to a full-factorial combination of three treatments: warming (ambient temperature or warming), flooding (no flooding, low-intensity flooding, or high-intensity flooding), and herbivory (no herbivory or herbivory), resulting in 12 treatment combinations. To exclude natural herbivory, we fenced both wetland areas during summer 2023. No flying herbivores (geese) or evidence of geese were ever observed in the enclosed areas. In three replicates of each treatment combination per wetland (36 mesocosms per wetland), we recorded air (+ 10 cm) and soil (− 5 cm) temperature every 90 min from early June to mid-August using dataloggers (iButtons DS1921G/Z, Maxim Integrated, San Jose, California, USA) (Appendix 1: Figs. S3-S4).

Warming treatment

To increase temperature, we installed ITEX-style (Hollister et al. 2022), conical open-top chambers (OTCs; 85 cm base × 50 cm top × 35 cm height; Kalwall Corp., Manchester, New Hampshire, USA) from 09-Jun-2023 to 11-Aug-2023. During the summer, OTCs increased mean air temperature by ~ 1.2 °C and soil temperature by ~ 0.2 °C. The increase in air temperature aligns with projected summer warming in the Y-K Delta, beyond warming that has already occurred, for the period 2040–2049 (SNAP 2020). At our site, OTCs had negligible effects on soil moisture, air relative humidity, and free-air CO_2_ concentration (Petit Bon et al. 2024).

Flooding treatment

To simulate flooding, we added high-tide water during three flooding events in mid-June, early July, and late July. The water source for the Lowland wetland was a small slough located approximately 50 m from the mesocosms, while in the Upland wetland we used a larger slough located about 30 m away.

We applied two flooding intensities to simulate different levels of high-tide inundation. Low-intensity flooding represented two consecutive days of overbank high-tide flooding per month, while high-intensity flooding represented six consecutive days of overbank high-tide flooding per month. These durations were based on projected future high-tide levels from relative sea-level rise forecasts for western Alaska (Sweet et al. 2022) and on the established relationship between tidal dynamics in our study area and the nearest tidal monitoring station in Nome (~ 330 km away; Terenzi et al. 2014). In the Lowland wetland, low-intensity flooding corresponds to conditions projected for the late 2020s, and high-intensity flooding to the mid-2030s. In the Upland wetland, low-intensity flooding is expected by the mid-2040s, and high-intensity flooding by the late 2050s (see Petit Bon et al. 2024 for further details).

On each flooding day that we measured, we added water until it reached the top of each mesocosm, on average 8.3 ± 2.6 SD L across all days and flooding events. We did not drain mesocosms because flooding treatments were meant to wet soils not only during flooding events but thereafter as would occur during a natural flooding event. Between flooding events, transpiration and evaporation reduced water levels to such an extent that we added approximately the same amount of water the first day of each flooding event, with rainfall likely accounting for differences. During flooding events, we added about 2–3 L to each mesocosm per day, signifying containers were more saturated during flooding events than between each flooding event, but there was significant loss each day to transpiration and evaporation even during flooding events.

We collected floodwater samples at the beginning and end of each consecutive flooding day and flooding event and calculated seasonal averages for total dissolved solids (1.08 ± 0.36 g L^−1^ in Lowland and 0.70 ± 0.64 g L^−1^ in Upland) and pH (6.8 ± 0.4 and 7.5 ± 0.3). These values were consistent with measurements from high-tide water samples collected daily throughout the summer (total dissolved solids: 0.84 ± 0.75 g L^−1^; pH: 7.1 ± 0.3) using a multiparameter water meter (HI98194, Hanna Instruments, Woonsocket, Rhode Island, USA).

Herbivory treatment

To apply the herbivory treatment, we simulated the three main components of goose herbivory: (1) grazing, (2) fecal deposition, and (3) trampling. The treatment was based on the intensity of herbivore pressure exerted by Branta nigricans (Pacific black brant), a common grazer in coastal wetlands. As low-lying coastal habitats become increasingly inundated, geese are expected to shift inland, applying similar grazing pressure to currently less-impacted areas such as our study wetlands. Therefore, our herbivory treatment represents current levels of herbivory in more heavily grazed coastal areas.

To simulate grazing, we clipped vegetation by removing leaves, buds, and fruits from both herbaceous and woody plants, as well as some herbaceous stems, until reaching the biomass target. Because goose diets can include substantial non-graminoid tissues, such as leaf tissue of other plants and berries from shrubs (Sedinger and Raveling 1984; Cadieux et al. 2005), we included these parts in our clipping protocol. Prior to clipping, we quantified the percent cover of each plant functional group (graminoids, deciduous dwarf-shrubs, evergreen dwarf-shrubs, and forbs) in each mesocosm using the point-intercept method (Bråthen and Hagberg 2004) with a 7 × 7 crossed grid (49 points). To reflect that more abundant species are likely grazed more frequently, we clipped each functional group in proportion to its relative cover (cf. Choi et al. 2022), thereby applying uniform grazing pressure relative to availability. This also allowed for meaningful comparison of plant trait responses among functional groups. We removed 26.4 g m^−2^ dry weight of aboveground biomass per mesocosm, corresponding to ~ 10% of peak-season biomass (Petit Bon et al. 2024), consistent with mean seasonal offtake by geese (Person et al. 1998). Clipping occurred twice during the season: once in late June and again in mid-to-late July. Biomass removal was greater in the second event (17.6 g m^−2^ vs. 8.8 g m^−2^) to reflect increased post-hatch grazing pressure, when goslings grow rapidly and adults molt in preparation for fall migration.

To simulate fecal deposition, we added 27 goose feces m^−2^ in mid-June, following the average natural deposition rate for B. nigricans throughout the summer (Beard et al. 2023). We collected fresh feces from the study area and distributed them evenly within each mesocosm. Feces were applied once early in the season to allow for prolonged nutrient effects and to ensure presence during trampling.

To simulate trampling, we applied pressure equivalent to the average body mass of B. nigricans (~ 1.2 kg), which exerts approximately 15.5 g cm^−2^ per step. We applied this load twice: after fecal addition and again following the second grazing event.

Focal species and functional trait sampling

We assessed aboveground plant trait responses to warming, flooding, and herbivory for the two dominant species in each wetland: C. rariflora and S. fuscescens in the Lowland wetland, and C. lyngbyei and S. fuscescens in the Upland wetland. This sampling design allowed us to compare (1) the magnitude of trait responses between species from different functional groups (sedges and deciduous dwarf-shrubs) within a given wetland, which are known to differ inherently in trait values (Chapin et al. 1996; Thomas et al. 2018), and (2) the response of a single species (S. fuscescens) growing under different biotic and abiotic conditions in the two wetlands.

Trait sampling was conducted between 31-Jul-2023 and 5-Aug-2023. From each mesocosm, we collected focal traits from two fully mature individuals of each focal species, selected at random, except for C. rariflora (Lowland wetland), for which we sampled three individuals. In total, we sampled 616 individuals. From each Carex individual, two leaves were collected; from each S. fuscescens individual, three leaves were collected. Leaves were pooled per individual for trait measurements. For S. fuscescens, leaves were selected from sun-exposed positions at or above the graminoid canopy, following standardized protocols (Pérez-Harguindeguy et al. 2013).

To isolate the direct effects of herbivory, we only sampled individuals that exhibited evidence of grazing in herbivory-treated mesocosms. For Carex spp., this meant leaves visibly clipped. For S. fuscescens, we sampled intact leaves from individuals that had other leaves removed during the herbivory treatment. This distinction reflects typical goose foraging behavior: graminoids are grazed at the leaf tip, leaving a portion intact, whereas shrubs with smaller, non-elongated leaves are more likely to experience full leaf removal.

We focused on four widely used aboveground plant traits, representing the two major axes of plant functional trait variation (Díaz et al. 2016). We measured two size-related traits, vegetative height (cm) and leaf area (mm^2^), and two resource economics traits, specific leaf area (SLA; mm^2^ g^−1^) and leaf dry matter content (LDMC; g g^−1^). While we acknowledge that leaf thickness, which can be more plastic than LDMC and vary independently, may provide complementary information on leaf structure, we focused on LDMC because it is related to thickness (Blumenthal et al. 2020) and the SLA-LDMC continuum is often used to capture the main axis of leaf economic strategies. Plant height was measured in the field on undisturbed individuals as the vertical distance from the ground to the highest vegetative tissue. After harvest, leaves were kept moist in zip-lock bags with dampened paper towels and refrigerated for a maximum of 48 h before processing. We first recorded fresh mass, then scanned leaves (Perfection V600 Photo Scanner, Epson America Inc., Los Alamitos, California, USA) and calculated leaf area using ImageJ (Schneider et al. 2012). Leaves were then oven-dried at 60 °C to constant mass and weighed. SLA was calculated by dividing leaf area by dry mass, and LDMC by dividing fresh mass by dry mass. Although LDMC is defined using water-saturated fresh mass (Pérez-Harguindeguy et al. 2013), leaf rehydration was not feasible under our field conditions at this remote sampling location. Yet, our approach ensures consistent comparisons within our experiment and aligns with procedures commonly used at remote tundra sites (e.g., Jónsdóttir et al. 2022; Kemppinen & Niittynen 2022).

Statistical analysis

We conducted all analyses in R v. 4.4.0 (https://www.r-project.org/). For each species, we averaged trait values across individuals within each mesocosm prior to analysis. We used the ‘vegan’ package to conduct permutational multivariate analysis of variance (PERMANOVA; Oksanen et al. 2020), the ‘prcomp’ function in base R for principal component analysis (PCA), and the ‘glmmTMB’ package (Brooks et al. 2017) to fit generalized linear mixed-effects models (GLMMs). We used the ‘car’ package (Fox & Weisberg 2018) for ANOVA and the ‘emmeans’ package (Lenth 2021) to calculate marginal means and perform pairwise comparisons. Model fit was assessed using the ‘DHARMa’ package (Hartig 2024), and lognormal transformations were applied when necessary to improve residual normality and homoscedasticity. We used the ‘ggplot2’ package (Wickham 2016) and the ‘factoextra’ package (Kassambara 2020) for data visualization.

To explore inherent differences in trait values among species of different plant functional groups and potential intraspecific variation in S. fuscescens across wetlands, we first tested for interspecific trait differences, treating S. fuscescens in the two wetlands as distinct entities. To assess multivariate trait differences, we performed a PERMANOVA using species as a fixed-effect on the centered and scaled Euclidean distance matrix of the four traits (height, leaf area, SLA, and LDMC). Row of mesocosms was included as a random-effect. Results were visualized using PCA. Complementarily, we tested each trait individually using GLMMs with species as a fixed-effect and row as a random-effect. We assessed species differences using Type III Wald chi-square tests via ANOVA and performed pairwise comparisons using Tukey-adjusted contrasts.

To evaluate treatment effects on multivariate trait expression, we conducted a separate PERMANOVA for each species and wetland type, using the same approach described above. Fixed-effects included warming, flooding, herbivory, and all two- and three-way interactions; row was included as a random-effect. Because none of the interactions were significant (P > 0.05), we report results from models including only the main effects. Multivariate trait responses were visualized via PCA, and loadings and scores of individual treatment effects were displayed separately for interpretability.

To examine treatment effects on individual traits, we used a separate GLMM for each trait, species, and wetland. Fixed-effects included warming, flooding, herbivory, and their interactions; row was included as a random-effect. When highest-order interactions were not significant, we sequentially removed the least significant term until either all retained interactions were significant or only main effects remained. We assessed treatment differences using Type III Wald chi-square tests via ANOVA and performed pairwise comparisons. For flooding (a three-level factor), we conducted pairwise contrasts between pre-specified levels (low-intensity vs. no flooding and high-intensity vs. no flooding); as only these planned comparisons were tested, no correction for multiple comparisons was applied (Gelman and Hill 2007). We present z-values and p-values from the conditional model to show all flooding levels.

Results

Interspecific functional trait differences

Plant species differed significantly in their multivariate trait expression (PERMANOVA: F_3,216_ = 60.2, P < 0.0001; Fig. 1, Appendix 1: Table S1). The PCA indicated that the sedge C. lyngbyei was associated with the greatest vegetative height and leaf area, C. rariflora with the highest LDMC and lowest SLA, and the deciduous dwarf-shrub S. fuscescens with generally lower height and leaf area, and higher SLA and lower LDMC compared to the sedge C. rariflora.Fig. 1. Plant species in the multivariate functional trait space. PCA biplot of plant traits: height, leaf area (LA), specific leaf area (SLA), and leaf dry matter content (LDMC), shown as black arrows. Each point represents the PCA scores of one species in a mesocosm and is colored by species identity. S. fuscescens is shown separately for the two wetlands. Ellipses indicate 95% confidence intervals around species centroids. Percent variance explained by the first two principal components is shown on the axes

Accordingly, species also differed significantly in individual trait values (Appendix 1: Table S2). The three species differed in height (P < 0.001) in the order C. lyngbyei > C. rariflora > S. fuscescens, with no height difference between S. fuscescens individuals in the two wetlands. All pairwise differences in leaf area were significant among species, including between S. fuscescens in the Lowland and Upland wetlands (P < 0.05), with C. lyngbyei having the largest leaf area, followed by C. rariflora, S. fuscescens in the Upland wetland, and S. fuscescens in the Lowland wetland. C. rariflora had lower SLA than all other species (P < 0.001), while C. lyngbyei and S. fuscescens (in both wetlands) did not differ in SLA. LDMC was higher in C. rariflora than in C. lyngbyei and S. fuscescens in the Lowland wetland (P < 0.001), but not different from S. fuscescens in the Upland wetland. S. fuscescens in the Upland wetland had greater LDMC than in the Lowland wetland, while C. lyngbyei did not differ from either.

Treatment effects on plant functional traits

The multivariate trait space of C. rariflora was not affected by warming, but was affected by flooding (PERMANOVA: F_2,51_ = 2.5, P = 0.07) and strongly by herbivory (F_1,51_ = 56.6, P < 0.001; Fig. 2A, Appendix 1: Table S3). However, at the single-trait level (Fig. 3A, Appendix 1: Table S4), warming increased leaf area by 15% (z = 2.1, P < 0.05). Low- and high-intensity flooding increased SLA by 8% (z = 2.8, P < 0.01) and 6% (z = 2.1, P < 0.05), respectively. Low-intensity flooding also decreased LDMC by 6% (z = − 1.9, P = 0.06). Herbivory reduced height by 35% (z = − 7.7, P < 0.001), leaf area by 50% (z = − 10.1, P < 0.001), and LDMC by 12% (z = − 7.3, P < 0.001), while increasing SLA by 14% (z = 5.6, P < 0.001).Fig. 2. Effects of treatments on multivariate functional trait expression. PCA biplots and PERMANOVA results of plant traits under warming, flooding (low- and high-intensity), and herbivory. Panels show results for: (A) C. rariflora, (B) S. fuscescens in the Lowland wetland, (C) C. lyngbyei, and (D) S. fuscescens in the Upland wetland. While treatment-specific biplots are displayed separately for interpretability, PCA and PERMANOVA were conducted once per species. Ellipses indicate 95% confidence intervals around treatment centroids. The four focal traits are shown as black arrows in the left column, which depicts the general trait structure for each species. Biplots are bolded when treatment effects were marginally significant (P < 0.10) and marked with asterisks when significant (P < 0.05). Percent variance explained by the first two principal components is shown on the axes in the left columnFig. 3Main effects of treatments on individual functional traits. GLMM predictions (± 95% confidence intervals) of plant traits under warming, flooding, and herbivory. Panels show modeled means for each species: (A) C. rariflora, (B) S. fuscescens in the Lowland wetland, (C) C. lyngbyei, and (D) S. fuscescens in the Upland wetland. Treatment levels include ambient temperature (A), warming (W), no flooding (nF), low-intensity flooding (LF), high-intensity flooding (HF), no herbivory (nH), and herbivory (H). Means and confidence intervals were back-transformed to the response scale for visualization. Bolded values indicate marginally significant effects (P < 0.10); asterisks denote significance at P < 0.05. Note that y-axis scales differ across species. For S. fuscescens in the Upland wetland, the interaction between flooding and warming on vegetative height is not visualized here; see Fig. 4 for details

For S. fuscescens in the Lowland wetland, the multivariate trait space was affected by warming (PERMANOVA: F_1,51_ = 2.3, P = 0.09) and strongly affected by herbivory (F_1,51_ = 24.6, P < 0.001), but not by flooding (Fig. 2B, Appendix 1: Table S3). At the single-trait level (Fig. 3B, Appendix 1: Table S4), warming increased height by 15% (z = 1.9, P = 0.06) and leaf area by 19% (z = 2.2, P < 0.05). High-intensity flooding decreased height by 15% (z = − 1.7, P < 0.10) and leaf area by 16% (z = − 1.7, P = 0.09). Herbivory decreased height by 17% (z = − 2.4, P < 0.05), leaf area by 36% (z = − 5.5, P < 0.001), and LDMC by 20% (z = − 5.3, P < 0.001), and increased SLA by 40% (z = 5.8, P < 0.001).

The multivariate trait space of C. lyngbyei was affected only by herbivory (PERMANOVA: F_1,45_ = 36.1, P < 0.001), with no effects of warming or flooding (Fig. 2C, Appendix 1: Table S3). At the single-trait level (Fig. 3C, Appendix 1: Table S4), herbivory decreased height by 30% (z = − 4.6, P < 0.001), leaf area by 44% (z = − 5.0, P < 0.001), and LDMC by 23% (z = − 7.3, P < 0.001), and increased SLA by 27% (z = 6.7, P < 0.001). Warming and flooding had no effects on any individual trait.

For S. fuscescens in the Upland wetland, the multivariate trait space was affected by flooding (PERMANOVA: F_2,53_ = 2.6, P = 0.03) and herbivory (F_1,53_ = 19.7, P < 0.001), but not by warming (Fig. 2D, Appendix 1: Table S3). However, at the single-trait level (Fig. 3D, Appendix 1: Table S4), warming decreased LDMC by 16% (z = − 1.9, P = 0.06). High-intensity flooding decreased SLA by 9% (z = − 1.8, P = 0.07) and increased LDMC by 11% (z = 2.4, P < 0.05). Herbivory decreased leaf area by 33% (z = − 4.7, P < 0.001) and LDMC by 16% (z = − 6.2, P < 0.001), and increased SLA by 25% (z = 5.2, P < 0.001). There was a significant warming × flooding interaction for plant height: flooding reduced height under ambient conditions, while warming increased height, but only under low-intensity flooding (z = 2.4, P < 0.05; Fig. 4).Fig. 4. Two-way interaction between flooding and warming on vegetative height of S. fuscescens in the Upland wetland. GLMM predictions (± 95% confidence intervals) under combinations of flooding and temperature treatments. Treatment levels include ambient temperature (A), warming (W), no flooding (nF), low-intensity flooding (LF), and high-intensity flooding (HF). Bolded values indicate marginally significant effects (P < 0.10); asterisks denote significance at P < 0.05

Discussion

Using a mesocosm experiment in two high-latitude coastal wetlands of the Y-K Delta, differing in landscape position and associated biotic and abiotic characteristics, we assessed the trait plasticity of three dominant plant species in response to multiple global change drivers. A single summer of simulated warming, flooding, and herbivory was sufficient to induce rapid changes in key aboveground functional traits, indicating a pronounced short-term capacity for phenotypic adjustment. These drivers acted largely additively, with limited evidence for interactive effects. However, the overall responsiveness and magnitude of trait plasticity varied across drivers and species. Herbivory consistently exerted the strongest influence on both size-related and leaf economics traits across all species, while responses to flooding and warming were more variable and dependent on species and wetland communities. This combination of consistent and context-dependent responses underscores the complexity of trait adjustment to multiple environmental changes and provides an important step toward predicting how high-latitude coastal wetlands may respond functionally to a changing climate.

Our finding that warming increased vegetative height (+ 15%) and leaf area (+ 15–19%) of C. rariflora and S. fuscescens in the Lowland wetland supports our prediction that higher temperatures promote increases in size-related trait values, indicative of greater biomass. This aligns with warming experiments showing positive effects on plant size in high-latitude ecosystems (Hudson et al. 2011; Baruah et al. 2017; Wei et al. 2023) and globally (Dobson and Zarnetske 2025), as well as with observational studies across latitudinal, elevational, and micro-topographical temperature gradients (Betway et al. 2021; Kemppinen and Niittynen 2022; Simin et al. 2022). However, no warming effects were detected in C. lyngbyei or S. fuscescens in the Upland wetland. Betway et al. (2021) similarly observed that while most tundra species increased in size across a 7 °C latitudinal gradient, others remained relatively stable. Stability may reflect inherently low trait plasticity but could also indicate “diminishing returns” on investment in size, whereby further allocation to structural growth yields limited functional benefits and constrains responsiveness to warming (Niklas et al. 2007). Such constraints may also help explain the limited response of C. lyngbyei to flooding. Similarly, Baruah et al. (2017) found that the deciduous dwarf-shrub S. herbacea increased in height in dry and tussock tundra but not in wet meadows when exposed to higher temperatures, suggesting that wetland context may have constrained S. fuscescens responses in the Upland wetland. In our study, this could be partly explained by competition or shading from C. lyngbyei, which is five times taller and produces ten times the leaf area of S. fuscescens. Despite predictions that warming would promote more acquisitive trait values, especially in the shrub, leaf economics traits showed minimal responses, except for a slight decrease in LDMC in S. fuscescens in the Upland wetland. While this weak response contrasts our expectations, it fits within broader mixed evidence: although increased acquisitive traits are the most reported response to warming across the Arctic (Bjorkman et al. 2018; Jessen et al. 2020; Wei et al. 2023), some studies have found shifts toward more conservative trait values (Hudson et al. 2011). Our findings suggest that leaf economic traits may be less responsive to short-term warming than size-related traits and are strongly shaped by species- and context-dependent plasticity (Hudson et al. 2011; Kemppinen and Niittynen 2022).

In line with our hypothesis, high-intensity flooding reduced vegetative height and leaf area of S. fuscescens in the Lowland wetland by 15% and 16%, respectively, consistent with negative impacts of inundation on shrub growth (Person and Ruess 2003; Churchill et al. 2015) and with broader patterns of height decline along flooding gradients (Baastrup-Spohr et al. 2015; Fu et al. 2015; but see Purcell et al. 2019). In the Upland wetland, both flooding intensities decreased S. fuscescens height under ambient but not warmed conditions, suggesting that warming may buffer flooding-induced growth reductions in certain contexts. Also, in line with our hypothesis, high-intensity flooding altered resource economics traits, increasing LDMC (+ 11%) and decreasing SLA (− 9%) in S. fuscescens in the Upland wetland, indicating a shift toward more conservative trait values under stress. Conversely, C. rariflora in the Lowland wetland exhibited a more acquisitive response, with increased SLA (+ 6%) and decreased LDMC (− 6%), the latter occurring only under low-intensity flooding. These contrasting species-specific responses mirror findings from lower-latitude systems, where flooding selects conservative or acquisitive strategies depending on species and context (Violle et al. 2011; Purcell et al. 2019), with plasticity further contributing to the dominant trait expression (Jung et al. 2010). In our system, the conservative response of the slower-growing dwarf-shrub contrasted with the more acquisitive response of the faster-growing sedge, which is better adapted to inundation. These findings support the potential role of flooding in reinforcing graminoid dominance in these high-latitude wetlands (Petit Bon et al. 2024). Notably, flooding effects on resource-use traits were evident in C. rariflora but absent in C. lyngbyei, and in S. fuscescens only in the Upland wetland, highlighting both species-specific responses, even among congeners, and the fundamental role of environmental context in shaping trait plasticity under new flooding regimes.

Herbivory, simulated as goose grazing, feces addition, and trampling, triggered a coordinated trait response across all three species, consistently affecting all four functional traits and emerging as the only driver to elicit a response in the otherwise unresponsive sedge C. lyngbyei. In line with our hypothesis, herbivory reduced vegetative height (17–35%) and leaf area (33–50%), aligning with known negative effects of grazing on plant size (Díaz et al. 2001; Barbero-Palacios et al. 2024). In sedges, these reductions likely resulted from direct defoliation, which in our study overcame any new growth typically observed following herbivory (Cargill and Jefferies 1984; Beard et al. 2019). In the deciduous dwarf-shrub S. fuscescens, the decline in size-related traits may instead reflect an avoidance response (Díaz et al. 2007; Jessen et al. 2020), and trampling, which could contribute to reduced height and constrained leaf development, as shrubs are particularly sensitive to physical disturbance (Egelkraut et al. 2020; Tuomi et al. 2021). Grazing may have also kept S. fuscescens in earlier vegetative stages. The observed increases in SLA (14–40%) and decreases in LDMC (12–23%), suggest that plants in these wetlands generally resist herbivory by investing in rapid regrowth (Díaz et al. 2001; Niu et al. 2016), consistent with our hypothesis. This also indicates for economic traits that S. fuscescens did not adopt a conservative strategy for herbivore avoidance. Nutrient input from feces may have further enhanced SLA, in line with fertilization experiments (Iturrate-Garcia et al. 2020; Jessen et al. 2020). Because high SLA is associated with rapid growth and high photosynthetic rates (Wright et al. 2004; Reich 2014), our results suggest that dominant species in high-latitude coastal wetlands maintain trait values that sustain performance and competitive ability under herbivory.

Warming, flooding, and herbivory primarily affected focal plant functional traits through additive rather than interactive effects. However, species varied in their sensitivity to these drivers. C. rariflora was highly responsive: flooding and herbivory together shifted this sedge toward more acquisitive trait values, with both drivers increasing SLA and decreasing LDMC, resulting in one of the stronger changes in leaf economics among all species. In contrast, warming and herbivory had opposing effects on size-related traits in both C. rariflora and S. fuscescens in the Lowland wetland, indicating that higher temperatures may partially offset the strong reduction in plant size caused by herbivory. These observations highlight the need to consider the additive effects of multiple global change drivers, which may reinforce or offset each other, in shaping the magnitude and direction of immediate trait change. Thus, even in the absence of interactions, each driver contributed significantly to trait variation, underscoring the need for multifactorial approaches to better predict the functional consequences of climate change in Arctic coastal ecosystems.

In summary, our one-year study demonstrates that warming, flooding, and herbivory – three major drivers of environmental change in high-latitude coastal wetlands of the Y-K Delta and beyond – can elicit immediate trait responses in dominant plant species. Herbivory reduced plant size and promoted acquisitive trait values, even in the deciduous dwarf-shrub S. fuscescens, which is thought to be less herbivore-tolerant than the two sedges (C. rariflora and C. lyngbyei). This pattern for the deciduous shrub may reflect the tendency of high-latitude deciduous plants to rapidly maximize carbon gain (Johnson and Tieszen 1976); shifting to a more acquisitive strategy under herbivory could indicate their ability to further accelerate growth when stressed. Flooding emerged as a stronger control than warming, reducing the size of S. fuscescens and shifting economics trait values in line with species-specific flood tolerance: toward conservative values in the less flood-tolerant S. fuscescens and acquisitive values in the more flood-tolerant C. rariflora. However, C. lyngbyei did not respond as C. rariflora did, and economics trait shifts in S. fuscescens were detected only in the Upland, not the Lowland, wetland. These results underscore the importance of understanding both intra- and interspecific trait responses across varying biotic and abiotic contexts to better predict future shifts in plant functional composition. Phenotypic adjustments at these levels may mediate plant community responses and ecosystem functioning even without species turnover (Siefert et al. 2015; Henn et al. 2018; Westerband et al. 2021; Jónsdóttir et al. 2022). While we focused on commonly used aboveground traits central to plant performance and biogeochemical processes (Bjorkman et al. 2018), belowground trait responses particularly to flooding, remain largely unknown (Purcell et al. 2019; Freschet et al. 2021; Pan et al. 2022). Expanding our work to include belowground trait dimensions will be essential for improving predictions of plant community and ecosystem responses to multiple, co-occurring global change drivers in coastal high-latitude wetlands.

Supplementary Information

Below is the link to the electronic supplementary material.Supplementary file1 (DOCX 8361 KB)