Aim: This study aims to evaluate whether the trait-environment relationships known for terrestrial plant height extend to freshwater macrophytes, and to identify the primary environmental factors influencing global height patterns in emergent and fully submerged aquatic plants.

Location: Global.

Time Period: Species occurrence records and height measurements compiled from 1931 to 2022.

Major Taxa Studied: Freshwater plant species.

Methods: We compiled a global dataset of maximum plant height records for 1,735 aquatic plant species, categorized by life form (partially emergent vs. submerged). Using generalized additive models, we tested how plant height varied along latitudinal gradients and how it was related to environmental predictors including temperature, relative inorganic carbon supply, and water depth (indicated by habitat availability). We explicitly analyzed whether these relationships differed between the two life forms, and validated these patterns using site-scale species and environmental data from northern temperate lakes and streams.

Results: Species with emergent growth increase in height with warmer temperature. On the contrary, fully submerged species exhibit increased height with higher inorganic carbon availability and colder climates. At the local scale, submerged plant height was positively associated with relative inorganic carbon supply in lakes and streams of the northern temperate zone.

Main Conclusions: The relationships between environmental variables and plant height differ between partially emergent to completely submerged life forms, consistent with realm-specific patterns in trait-environment correlations. These differing patterns highlight that trait-environment coupling in freshwater systems may not parallel those in terrestrial ecosystems and suggest that life form mediates species’ responses to environmental variation.

2.1 | Data compilation

To study global patterns of aquatic plant height, we searched for height data for all known aquatic plants occurring in freshwater (3696 species; summarized by Murphy et al. (2019)). We collated height information for roughly half of the world’s aquatic plant species, i.e., 1,735 species (1,206 emergent plants, 411 submerged plants, 128 floating plants). The data comprised 6,186 records of height data and life form information for aquatic plants from 147 references and 120 online databases (including 810 records from the TRY database (Kattge et al. 2020), Supporting information Appendix 1) accessed from 22^nd Dec 2021 to 20^th Jan 2022. For each species, our dataset consists of an average of 3.57 (1 - 163) records. Where multiple records were available, we selected the maximum reported height to represent the species’ potential height under favorable conditions. Although we aimed to use vegetative height, this distinction was not consistently provided in some original sources (e.g. TRY). Therefore, our maximum height values may sometimes include reproductive structures. We acknowledge this as a limitation of our dataset, however, we believe this does not introduce systematic bias across life forms, and the broad patterns we report remain robust. For species with unknown life form (54 species), we used visual inspection of species photographs for categorizations (submerged vs. emergent life form) based on the location and types of leaves of the species. Species were classified as emergent when more than half of their leaves were positioned above the water surface, and as submerged when more than half of their leaves were entirely below the water surface. For a subset of the submerged species with recorded plant height, their photosynthetic usage of inorganic carbon was known as either strictly CO₂ (40 species) or both bicarbonate and CO₂ (39 species, Data S4; Iversen et al. 2019), thus allowing for a division into bicarbonate-users and obligate CO₂-users. This set of species was used to test whether height differed between the two inorganic carbon strategies.

We combined the collected height data with the presence/absence information for each species in 10° × 10° (latitude × longitude) grid cells at global scale extracted from Murphy et al. (2019), who compiled a global, grid‐based distribution dataset for aquatic macrophytes by synthesizing information from regional floras, monographs, published literature, and databases. By using the maximum height value recorded for each species, we calculated the average maximum height of the plant assemblage in each grid cell spanning the entire globe (hereafter abbreviated as grid cell plant height, which effectively provides an estimate of the average plant maximum height of the species with ranges overlapping a given grid cell). The grid cell plant height was calculated using the arithmetic mean of the height of the species in each assemblage. By adopting the arithmetic mean of species’ maximum heights as our grid cell-level metric, we capture the community’s average potential for growth and competitive performance, serving as a useful proxy in assessing plant ecological strategy (Westoby et al. 2002; Osada et al. 2014; de Vries et al. 2024). This resulted in a dataset of 237 grid cells spanning all vegetated land areas of the earth (Data S1). To quantify the patterns for the different life forms of aquatic plants, we separated the dataset into emergent species (occurring in 237 grid cells), submerged species (n=233 grid cells; Data S1). Floating plants were excluded in the analysis because their height to a large degree is structured by local variations in water depth (floating rooted), or they are often very small in size (free-floating: e.g., Wolffia arrhiza: 0.1 cm).

We conducted an additional analysis using available intraspecific height data to verify whether maximum height values reliably reflect a species overall trait profile and to identify any potential biases from using extreme values. Specifically, we applied linear regression on log-transformed species height data to examine the relationship between maximum and average heights. The results, showed a strong correlation between the true maximum recorded heights and the simulated species height values derived from repeated random draws (Pearson R = 0.98, p < 0.001). The mean absolute error of 0.56 log units between log-transformed simulated and true maximum heights further supported the robustness of our maximum height metric (Fig. S1).

2.2 | Environmental conditions

Climatic variables related to temperature and precipitation (see Appendix 2) were extracted from the WorldClim database (http://worldclim.org/version2), using data at 10-arc min resolution (Hijmans et al. 2005). We condensed temperature and precipitation variables using a PCA for each set of variables and extracted the 1^st-axis (PC1) as a proxy of temperature and precipitation, respectively. The 1-axis from the temperature PCA (PC1_temp) correlated strongly with most variables related to average temperature conditions such as the mean temperature in the warmest quarter (R = 0.94, Fig S2), and likewise the precipitation axis (PC1_precip) was strongly related to annual precipitation (R = 0.99, Fig S2).

As an indicator of the availability of shallow and deeper aquatic habitats, we collected data on the amount of seasonal (shallow) and permanent (deeper) water cover per grid cell from Pekel et al. (2016). We used the fraction of the total water covered area comprised by seasonal water (seasonal water ratio) to estimate dominance of either shallow or deeper freshwater habitats within each grid cell.

As the photosynthesis and growth of submerged plants can be limited by inorganic carbon supply due to the low diffusivity of CO₂ in water (Maberly & Gontero 2018), we also collected the inorganic carbon concentration for each grid cell. For CO₂, we estimated the average concentration in the grid cells using the global dataset (pCO₂; µatm from lakes and streams separately) from Raymond et al. (2013). For each grid cell, we calculated the CO₂ values by extracting mapped pCO₂ from lakes and streams and computed a weighted mean, where each pCO₂ value was weighted by the relative area of its corresponding waterbody within each COSCAT regions (Meybeck et al. 2006; Raymond et al. 2013) using ArcGIS Pro (version 2.7.0). For bicarbonate, we used the global bicarbonate dataset (HCO₃^-; mM) provided by Iversen et al. (2019) and analyzed it using the mean value for each grid cell.

To understand the role of the supply of inorganic carbon (CO₂ and HCO₃^-) to submerged plants, we estimated a supply index using a framework equivalent for estimating oxygen limitation for aquatic ectotherms (oxygen supply index, Verberk et al. (2011)). Modified from the oxygen supply index equations (see Appendix 3), we calculated the bicarbonate supply index (HCO₃SI; mol* m⁻¹ * s⁻¹) and the CO₂ supply index (CO₂SI; mol* m⁻¹ * s⁻¹) incorporating both the concentration of each compound and the temperature-dependent diffusivity. The sum of the two indices estimates the total supply of inorganic carbon (inorganic carbon supply index; ICSI; mol* m⁻¹ * s⁻¹). As the amount of sequestered carbon used for basic metabolism is expected to be dependent on temperature (Gillooly et al. 2001), we calculated the relative supply of inorganic carbon as the ICSI divided by the temperature-dependent carbon demand as suggested by Verberk et al. (2011) (relative inorganic carbon supply index, rICSI; mol* mol* m⁻¹ * s⁻¹) :

Where Q₁₀ is the temperature coefficient (using a standard value of 2.0, (Rasmusson et al. 2019)) and ΔT is the difference between the grid cell temperature (temperature in the warmest quarter) and the average temperature in the entire dataset (℃). Environmental variables are matched with grid-cell plant height data using grid-cell information and compiled in Data S2.

As a supplement to the pattern of grid cell plant height of submerged species at the global scale, we also analyzed a finer-scale dataset with site-specific records of plant species and measurements of CO₂ and bicarbonate from 963 sites across the Northern Hemisphere (lakes and streams from Europe and North America (Iversen et al. 2019)). Species occurrences were compiled from regional aquatic vegetation monitoring programs and previously published floristic surveys. Based on the extracted dataset, species richness per site ranged from 5 to 52 species, with a mean of 10 species per site.

We analyzed the relationships between the average maximum plant height (abbreviated as site plant height) and two environmental variables calculated at site level; relative inorganic carbon supply and a temperature variable (PC1_{site temp}) constructed in the same way as for the global data. After removing sites with incomplete data (e.g., unknown site names or missing species height records), the final dataset comprised 539 sites (Data S3).

2.3 | Statistical analysis

To account for phylogenetic constraints on plant height, we tested whether the maximum plant height contained a phylogenetic signal by estimating Blomberg’s K with the function ‘phylosignal’ in R package Picante (Blomberg et al. 2003). The phylogenetic relationships were obtained using phylomatic-awk-1.1.0 (Zanne et al. 2014). We extracted phylogenetic information for the selected species (1644 species matched) from the backbone tree in Zanne et al. (2014), which was constructed using multi-gene molecular and fossil data. Our analysis indicated that a very low phylogenetic signal in aquatic plant height (Blomberg’s K = 0.026), with a significant deviation from expectations under Brownian motion (p < 0.001). This suggests that plant height is not strongly structured by shared ancestry and may be influenced by other ecological or evolutionary factors (Blomberg et al. 2003). As a result, phylogeny was not taken into account in subsequent analyses.

To assess if there was a latitudinal pattern in aquatic plant height (including non-linear relationships), we first used Generalized Additive Models (GAM, mgcv package in R) to model the effect of latitude on grid cell plant height, with latitude as a smoothed term (Joswig et al. 2022). This model was compared to a constant null model using an F-test. Non-significance of the smoother term would indicate a random pattern between latitude and plant height. GAMs were selected because they allow flexible non-linear responses, which are expected along broad biogeographic gradients where biological patterns vary non-linearly with environmental and climatic drivers. Subsequently, we tested for potential divergent patterns in emergent and submerged plants by introducing plant life form and its interaction with latitude as independent variables in the model. To assess latitudinal patterns in plant height while accounting for spatial dependence, we fitted a generalized additive model (GAM) including a two-dimensional spatial smoother (s(Longitude, Latitude, bs = "sos")) together with an interaction between absolute latitude and life form. This approach controls for spatial autocorrelation in the data while allowing the latitudinal effect to differ between emergent and submerged species. To assess whether height variability differed between life forms at the spatial scale of analysis, we compared the variance in plant height between emergent and submerged species using Levene’s test for homogeneity of variance.

To examine whether emergent and submerged species differed in their responses to environmental gradients, we included interaction terms between life form and each environmental predictor in the GAM including the linear effect of seasonal water ratio, PC1_temp, PC1_precip and rICSI as explanatory variables. To account for variation in sampling robustness among grid cells, each grid cell was weighted by species richness in the model. We also verified that the results were robust to omitting the weighting scheme, with parameter estimates and significance remaining largely unchanged. To account for spatial autocorrelation, we included longitude and latitude in the model as a spatial spline. Due to a strong correlation between PC1_precip and rICSI (Pearson's correlation coefficient > 0.7), we chose to retain the environmental indicator rICSI, which we consider more ecologically more important, and removed PC1_precip from the model. rICSI was log-transformed to meet the assumption of a log linear effect on plant performance (Iversen et al. 2019) and all variables were subsequently normalized as Z-scores to ensure that parameter estimates were on the same scale to improve interpretation. The parameter estimate of this interaction term represents the differences in regression slopes between emergent and submerged species of a variable of interest. Significance was based on the 95% confidence interval of estimated parameter means.

Site plant height data (containing only submerged plants) was analyzed with a GAM model that included a spatial spline (latitude and longitude), using rICSI and temperature (PC1_{site temp}) as explanatory variables.

A Welch t-test was performed to test the height difference between the two inorganic carbon use strategies. All statistical analyses were performed in R software (version 4.1.3).