Introduction

The ongoing biodiversity loss due to human actions is an undeniable crisis (IPBES, 2019; IUCN, 2019). Land use changes, habitat degradation and introduction of invasive species to ecosystems are among the main causes of biodiversity loss (IPBES, 2019). One consequence of these anthropogenic changes is biotic homogenisation, a process where ecosystems lose their biological uniqueness, and similarity among communities increases (McKinney & Lockwood, 1999; Olden & Rooney, 2006). Contrary to this process, there is biotic differentiation, where similarity among communities decreases (Olden & Poff, 2003). Both processes are multidimensional, as they cover the loss or gain not only of taxonomic distinctiveness over time but also of functional, genetic (e.g. Olden et al., 2004) and phylogenetic (e.g. Harrison et al., 2018; Liang et al., 2019) distinctiveness. However, most studies have focused only on taxonomic homogenisation and differentiation, while functional and genetic aspects have received less attention (Olden et al., 2018; for definitions, see Table 1).

Table 1 Definitions of the general terms used throughout this paper
Full size table

The importance of functional traits and niche processes has been acknowledged for a relatively long time (e.g. McGill et al., 2006), but functional similarity and uniqueness have only recently received more attention (e.g. Gámez-Virués et al., 2015; Liang et al. 2019). Functional homogenisation, which refers to a loss of specialised species or entire functional groups (Olden et al., 2004), is directly linked to ecosystem functions and indirectly to ecosystem services through the loss of functional diversity (Clavel et al., 2011). The maintenance of multiple functions in an area thus needs multiple species with a variety of traits. It has been suggested that the decline of specialist species can cause global functional homogenisation (Clavel et al., 2011) and that biotic homogenisation is connected with ecosystem multifunctionality (e.g. Hautier et al., 2018). However, the role of beta diversity in ecosystem functioning is still not clear (Mori et al., 2018), and it is still unclear how functional homogenisation and differentiation are acting at different levels of biotic communities and under different levels of human impacts.

It has been argued that we are heading towards an era called ‘Homogocene’, as many studies have shown that human actions are causing increasing biotic homogenisation in different environments (Olden et al., 2018). As there are several studies showing that homogenisation is particularly strongly acting in aquatic environments (e.g. Rahel, 2002; Petsch, 2016), the term ‘Aquatic Homogocene’ has also been introduced. Studies concerning aquatic homogenisation have focused mainly on fish (e.g. Castaño-Sánchez et al., 2018; Richardson et al., 2018) and macroinvertebrates (e.g. Cook et al., 2018; Zhang et al., 2018b; dos Bertoncin et al., 2019), while other organism groups, such as macrophytes, have received less attention. Studies focusing on the biotic homogenisation of aquatic macrophyte communities have found contradicting results, showing either biotic homogenisation or differentiation (e.g. Lougheed et al., 2008; Johnson & Angeler, 2014; Baastrup-Spohr et al., 2017; Elo et al., 2018; Zhang et al., 2018a). These studies have used different approaches and have examined a variety of diversity measurements. Previous studies have been based on historical datasets and re-sampling the same sites (e.g. Baastrup-Spohr et al., 2017; Zhang et al., 2018a), one time point (e.g. Lougheed et al., 2008; Johnson & Angeler, 2014; Elo et al., 2018) or palaeoecological methods (e.g. Salgado et al., 2018).

Functional traits have been studied to some extent with aquatic macrophytes (e.g. Arthaud et al., 2012; Baattrup-Pedersen et al., 2016). However, these investigations have not been done extensively in a temporal framework. Functional homogenisation and differentiation of aquatic macrophyte communities have received even less attention than taxonomic ones. Only Zhang et al. (2018a) studied functional richness and compositional dissimilarities with aquatic macrophytes from before 1970s to after 2000s, and they found functional differentiation of macrophyte assemblages instead of homogenisation, even though their study area in the Yangtze River floodplain has faced strong human impacts for a long time. Furthermore, to our knowledge, there are no studies of functional changes in aquatic macrophyte communities at the species level.

Our aim was to focus on the issue of functional homogenisation and differentiation by examining if functional features of vascular aquatic macrophytes have changed or remained the same after a period of 70 years in boreal lakes (Fig. 1). To increase understanding of functional patterns and underlying mechanisms, we explored these patterns at the community and species levels. We used data of vascular aquatic macrophytes from 1947–1950 and 2017 based on surveys of the same set of boreal lakes. At the community level, we examined if there has been functional homogenisation or differentiation or no change in functional similarity. Also, at the community level, we examined if aquatic macrophyte communities show different spatial patterns in functional composition in relation to the environment between the two time periods. We also studied if changes in functional richness in time are linked to concomitant changes in the environment. To further understand the changes in community composition, we also ordered each species in functional space to see if species with certain sets of traits become more common or rare in the study period (representing the species level). As human impact has not been very strong in our study area, we predicted that there might not be drastic changes in functional features. However, we hypothesised that (1) different environmental variables explain functional community composition between the decades, (2) changes in functional richness are linked to changes in the environment across decades (e.g. Zhang et al., 2018a) and (3) declining species (i.e. ‘losers’) have different functional traits than increasing (i.e. ‘winners’) or stable species (e.g. Steffen et al., 2013).

Fig. 1
figure 1

Flow chart showing main questions, methods, results and conclusions. See text and Supplementary material Appendix 4 for detailed method and data descriptions

Full size image

Materials and methods

Study area

We studied 28 lakes located in the southern boreal zone (Ahti et al., 1968) in Finland, in the Kokemäenjoki drainage basin between the two large lakes Roine and Pyhäjärvi (Fig. 2; Supplementary material Appendix 1). Many of these lakes have brown, humic waters and are situated in small chains of lakes and streams, both features typical of boreal lakes. The elevation among the lakes studied varies from 77 to 131 metres above sea level. Nonetheless, changes in environmental conditions occur even along a relatively low elevation gradient in the study area, mainly due to glacial and postglacial processes. Due to these historical processes, fine-grained sediments with nutrients have been washed along the elevational gradient in the landscape, as is common in glaciated areas (Seppälä, 2005). Thus, naturally more eutrophic lakes are located at lower elevations and are mainly surrounded by agricultural land or settlements, while smaller and more oligotrophic lakes at higher elevations in the landscape are less affected by human activity and are surrounded by coniferous forest and peatlands. Key characteristics of these lakes are described in Table 2, and more information of the study area can be found in Lindholm et al. (2019).

Fig. 2
figure 2

Map showing the location of the study area and the lakes (N = 28) studied. Names of the lakes can be found from Supplementary material Appendix 1

Full size image
Table 2 Environmental variables, species richness values and descriptive statistic from the 1940s and 2017. Wilcoxon signed rank test results between environmental variables from 1940s and 2017 (N = 28)
Full size table

Species and trait data

We used aquatic macrophyte data from 28 lakes (Supplementary material Appendix 1) from two different time periods. The first macrophyte survey was conducted in 1947–1951 by Dr. U. Perttula (1954, unpublished), and the second survey was done using similar methods in the summer 2017 (July and August). Two different types of rakes and an aquascope were used in surveying aquatic macrophytes in the whole lake area. In this study, we concentrated on presence-absence data of species classified traditionally as aquatic vascular plants in Finland (Linkola, 1933) and, in addition, from the sedge genus Carex seven tall species that grow in the water. With historical datasets, presence-absence data are usually the most reliable source of information compared to coverage or other information representing abundance in macrophyte studies. We did not include taxa identified to genus level (e.g. Isoëtes sp.) or hybrids in the analysis. However, we combined some species to species complexes due to identification differences between the two decades (Supplementary material Appendix 2). Thus, in total, 65 vascular plant taxa were included (Supplementary material Appendix 2). Species richness values and descriptive statistic in the two decades can be found in Table 2.

To represent functional diversity of aquatic macrophytes, we used four functional traits: growth form, perennation, common method of propagation and potential size (Supplementary material Appendix 2; Lindholm et al. 2019). Growth form division consists of the following classes: ceratophyllid, elodeid, helophyte (incl. tall Carex), isoetid, lemnid and nymphaeid (Toivonen & Huttunen, 1995). Only the main growth form of species was considered. Perennation consisted of three ranked classes: (1) annual, (2) biennial/short lived perennial and (3) perennial. Perennation information was mainly collected from Willby et al.’s (2000) attribute-based data, where some species had an attribute present at two different categories. In such cases, the species obtained a value in between the ranked categories (i.e. 1.5 and 2.5), following Göthe et al. (2017). However, we decided to weight value 2, which indicates the presence of the attribute, at the expense of value 1, which indicates the occasional, but not general exhibition of the attribute. Perennation information was not available in this source for all species, and in such cases, data were complemented by information from other literature sources and databases (e.g. Ecological Database of the British Isles: http://ecoflora.org.uk/). Degree of vegetative propagation (Klotz et al., 2002) consists of five ranked classes: (1) by seed, (2) mostly by seed but also vegetatively, (3) by seed and vegetatively, (4) mostly vegetatively and also by seed and (5) vegetatively. Potential size information (cm) is a continuous trait from Hämet-Ahti et al. (1998) complemented from Mossberg and Stenberg (2012) for a few species. It represents the potential length of an individual omitting the root or rhizome length (Bornette et al., 1994; Doledec & Statzner, 1994).

Environmental variables

We used the following environmental variables: lake area (ha), elevation (m), maximum depth (m), pH, water transparency (m), land use variables and geographical coordinates of lakes’ centres, as we had only a limited number of environmental variables available from the 1940s (Table 2; Lindholm et al., 2019). Measurements of pH were done in the 1940s in summertime, while in 2017 during the fall overturn. Water transparency (m) was measured using a Secchi disk at the same time as the macrophyte sampling. Land use variables from 200 m buffer zones (Pedersen et al., 2006) were derived from the base maps for both decades (National land survey of Finland, 2017, 2018). We calculated three land use variables proportional to buffer area: agricultural area (i.e. field and pasture area), built area and amounts of ditches. These land use types have changed most over the past decades in the study area. We conducted the Wilcoxon signed rank test to detect differences in environmental variables between 1940s and 2017.

Data analysis

A schematic diagram showing the statistical methodology used can be found in Supplementary material Appendix 3. We did the following analyses to the time periods 1940s and 2017 using the vascular aquatic macrophyte presence-absence data and the functional trait data. All analyses were conducted in the R environment (R Core Team, 2017).

First, to explore functional homogenisation we calculated the amount of functional beta diversity (reflecting both species replacement and loss/gain) among all pairwise comparisons of lakes (Legendre, 2014) for the two time periods separately. We used the Gower distance (Gower, 1971) to calculate between-species distances based on the functional trait data using the function gowdis in the FD package (Laliberté & Legendre, 2010; Laliberté et al., 2014). Then, we used a hierarchical cluster analysis to this species-by-species matrix to produce a trait tree by using the function hclust in the stats package (R Core Team, 2017). Then, using the trait tree and macrophyte presence-absence data (site-by-species matrix), we calculated both the average and the variance of functional beta diversity based on Sørensen dissimilarity index (Legendre & Legendre, 2012) using the function beta.multi in the BAT package (Cardoso et al., 2015, 2018).

Second, we examined if aquatic macrophyte communities showed different spatial patterns among lakes in functional composition in relation to the environment between the 1940s and 2017. We did the following analyses separately for the time periods. We produced a dissimilarity matrix based on the same trait tree and the site-by-species matrix as mentioned above. We generated the functional dissimilarity matrix based on the Sørensen dissimilarity index (Legendre & Legendre, 2012) using the function beta in the BAT package (Cardoso et al., 2015, 2018). Then, we used distance-based redundancy analysis (db-RDA) to examine functional community-environment relationships (Legendre & Anderson, 1999). In this constrained ordination technique, the dissimilarity data are first ordinated using metric multidimensional scaling (principal coordinates analysis, PCoA), and the ordination results are analysed by RDA (Legendre & Anderson, 1999). This method works with any type of distance matrix as the response matrix (Legendre & Legendre, 2012). We selected predictor variables for the model using the forward selection procedure with two stopping rules (Blanchet et al., 2008) using the function ordiR2step in the vegan package (Oksanen, 2017). Then, we ran the db-RDA using the function capscale with the lingoes correction to avoid negative eigenvalues using the vegan package (Oksanen, 2017). Then we tested the model significance by 999 permutations using the function anova in the vegan package (Oksanen, 2017).

Third, we studied how the temporal change in functional richness (i.e. the proportion of the functional space filled by community) is linked to the changes in the environment. As we wanted to measure the total functional range covered by the community (Legras et al., 2018), we used functional richness index (FRic) of Villéger et al. (2008). We calculated FRic based on raw data using the function dbFD with square root correction in the FD package (Laliberté & Legendre, 2010; Laliberté et al., 2014) to each lake for the 1940s and 2017, respectively. Then, we calculated the change in FRic between the two time periods (from now on FRicc). We ran Spearman correlation tests between FRic and species richness in the 1940s and in 2017, and between the FRicc and the change in species richness between decades. We divided our environmental variables to ‘stable’ and ‘non-stable’ variables, as we wanted to find out if changes in FRicc are related to concomitant changes in the environment or environmental factors that are stable through decades. We calculated the change in non-stable environmental variables (pH, water transparency, agricultural area, built area and ditches) between the 1940s and 2017. For clarity, we refer to these changes as pHc, water transparencyc, agricultural areac, built areac and ditchesc. We used values from the 1940s with the more or less stable variables (east and north coordinates, elevation, lake area and depth). Next, we modelled the relationship between FRicc and all environmental variables using linear regression (LR). FRicc was log-transformed to normalise distribution. Collinearity was tested by variance inflation factor (VIF). The LR optimization was based on Akaike’s Information Criterion (AIC) (Burnham & Anderson, 2004).

In addition, we executed commonality analysis to decompose the linear regression R2 to unique and common components of explanatory variables. Unique components indicate how much variance is uniquely accounted for by a single variable, and common components how much variance is common to a variable set. We performed the commonality analysis by using the function commonalityCoefficients in the yhat package (Nimon et al., 2013). To further disentangle the relationship of functional richness to the environment, we examined FRic in the 1940s and 2017 separately. We forced into both LR models four predictor variables based on the previous LR model of FRicc, as the purpose was to compare the same predictor variables between the time periods.

Fourth, we disentangled the changes in macrophyte communities by observing each species in functional space, to see if species with certain sets of traits become more common or rare during the study period. Functional space is a multidimensional space, where the axes are combinations of functional traits and species are placed in the coordinates given by their functional traits (Carmona et al., 2016). We used the Gower distance (Gower, 1971) to calculate between species distances based on the trait data using the function gowdis in the FD package (Laliberté & Legendre, 2010; Laliberté et al., 2014). Then, we used ordination by non-metric multidimensional scaling (NMDS) using the function metaMDS in the vegan package (Oksanen, 2017). In the produced functional space, we plotted species based on NMDS with two dimensions (final stress level = 0.179). We represented the species with different symbols based on whether the species have (1) declined, (2) remained the same (stable) or (3) increased in distribution during the study period in the study area (i.e. differences in the number of localities between two different decades). A species was considered declining when its occurrence decreased by two or more lakes, and vice versa for increasing species. Stable species declined or increased in a maximum of one lake (i.e. from − 1 to 1) between the time periods. Finally, we used permutational multivariate ANOVA (PERMANOVA: Anderson, 2001) to test for significant differences between these three species groups by 999 permutations using the function adonis in the vegan package (Oksanen, 2017). We also associated NMDS species scores to each trait by conducting the Kruskal–Wallis rank sum test (Hollander & Wolfe, 1973) using the function kruskal.test in the stats package (R Core Team, 2017) and by plotting the relationships between the NMSD axes and the four traits.

Results

Changes in environmental variables

The Wilcoxon signed rank test results between the environmental variables from 1940s to 2017 can be found in Table 2. The changes in the non-stable environmental variables from the 1940s to 2017 can be found in Fig. 3. The changes in pH between time periods have been mixed, and the water transparency (m) in the lakes near the urban areas has increased. The lakes where water transparency has declined are located at higher elevations in the landscape, are less affected by human activity and are surrounded by coniferous forest and peatlands. Agricultural area showed a clear declining trend, while built area and ditches have increased in most lake areas (Fig. 3).

Fig. 3
figure 3

Maps showing the changes in the environmental variables and species richness from the 1940s to 2017 in the lakes (N = 28) studied. Elevation model from National Land Survey of Finland (2015) and water bodies from Finnish Environment Institute (2014)

Full size image

Functional dissimilarity

There were no clear changes between the 1940s and 2017 in functional beta diversity based on Sørensen dissimilarity index (Fig. 4). The average value of Sørensen dissimilarity index was 0.29 in the 1940s and 0.30 in 2017.

Fig. 4
figure 4

Variation in functional beta diversity based on Sørensen dissimilarity index in the 1940s and 2017. N = 28 lakes in both 1940s and 2017

Full size image

Changes in functional community composition

The final db-RDA detected three (P < 0.05) environmental variables related to the variation in functional community composition among the lakes both in the 1940s and 2017 (Fig. 5). In the 1940s, the variables that best accounted for the variation in functional community composition were elevation, depth and ditches. In 2017, pH, depth and area were the most important variables. The amount of explained variation (adj. R2) in both time periods was almost the same, 19.9% in the 1940s and 20.4% in 2017.

Fig. 5
figure 5

Plots of distance-based redundancy analysis of functional composition of aquatic macrophyte communities in the 1940s (A) and in 2017 (B). The significant environmental variables are shown as arrows in the plots. Names of the lakes can be found from Supplementary material Appendix 1

Full size image

Changes in functional richness

For the set of traits considered, dimensionality reduction occurred and only the first four PCoA axes were used to calculate FRic for the 1940s. As a result, the quality of the reduced-space representation of FRic was 0.69. For 2017, the first five PCoA axes were used and the quality of the reduced-space representation of FRic was 0.76. FRic in the 1940s was generally higher than in 2017 (Supplementary material Appendix 4, Fig. S1). The Spearman correlation between FRic and species richness in the 1940s was r = 0.92 (P < 0.001) and r = 0.91 (P < 0.001) in 2017. The correlation between the FRicc and the change in species richness between decades was 0.34 (P = 0.075) (Supplementary material Appendix 4, Fig. S2).

The final LR model included only one non-stable variable, agricultural areac, and three stable variables: lake area, elevation and north coordinate (Table 3). The LR model explained 67% of the variation in the change of the FRic (adj. R2 value) and the AIC value was 54.4. The assumptions of normal errors and independency were not statistically violated based on the exploration of residuals, and VIF for all explanatory variables was lower than three. The lakes at the lower elevation in the landscape have had more changes in functional richness between the time points than the lakes in the upper parts of the landscape (Supplementary material Appendix 4, Fig. S3). The size of the lake and the changes in functional richness showed a positive relationship, thus changes in the functional richness are higher in the larger lakes, even though the pattern was not as clear as with elevation. More information about the relationship between FRicc and the change in agricultural area can be found in the Supplementary material Appendix 4 (Figs. S4, S5).

Table 3 Summary of the results of linear regression (LR) analysis (response variable = change in functional richness index between the 1940s and 2017) and commonality analysis for each environmental variable based on linear regression model adjusted R2 values
Full size table

Based on both the unique and total effects of explanatory variables decomposed from linear regression (R2), elevation contributed the most to the whole LR model. However, the common contribution of elevation to the LR model was higher than the unique, indicating shared effects with the other variables. North coordinate contributed the second most to the whole model based on both the unique and total effects. The only non-stable variable included in the LR model contributed the least based on both the unique and total effects, and the common effect was slightly higher than the unique effect, indicating also shared effects with the other variables (Table 3).

The LR model in the 1940s explained 61% of the variation in the log FRic (adj. R2 value). In the 1940s, the only statistically significant variable was elevation. The LR model in 2017 explained 70% of the variation in the log FRic (adj. R2 value). In 2017, statistically significant variables were agricultural area, lake area and north coordinate (Supplementary material Appendix 4, Table S1).

Species in functional space

For each species in the functional space, there was no clear pattern of species with certain sets of traits becoming more common or rare during the study period (Fig. 6). In addition, based on PERMANOVA, there was no significant difference between the groups of declined, stable or increased species (F = 1.402, 62, P = 0.243). The four traits were significantly different from either one of the NMSD axes based on Kruskal–Wallis rank sum test (Supplementary material Appendix 5).

Fig. 6
figure 6

Non-metric multidimensional scaling (NMDS) plots of species distribution in functional space. In plot A, shown are species abbreviations (Supplementary material Appendix 2), and in plot B, the lines connect each species to its group centroid. Green = stable species, blue = increased species and grey = decreased species

Full size image

Discussion

In many geographical areas and organism groups, the functional homogenisation process has been noticed to occur, thereby seriously affecting ecosystem functioning (e.g. Bergeron et al., 2019). However, we did not find signs of functional homogenisation or, on the other hand, functional differentiation in our study area during a 70-year period for vascular aquatic macrophytes, as the multiple-site functional beta diversity remained almost unchanged. Our prediction that there might not be drastic changes in the functional features was supported, even though we detected some changes. Our hypothesis that different environmental variables would explain the functional community composition between the decades was partly supported. We found that changes have occurred in the functional community composition-environment relationship, as the aquatic macrophyte communities showed different spatial patterns in functional composition between the two time periods, and the main environmental drivers for these changes were partly different. Our second hypothesis that the changes in functional richness would be linked to changes in the environment across decades (e.g. Zhang et al., 2018a) was also partly supported, as we found that changes in functional richness in time were only partially linked to the concomitant changes in the environment. Stable factors, mainly those related to the position of the lakes in the landscape, were more important in explaining the changes in functional richness. Our third hypothesis that declining species (i.e. ‘losers’) would have different functional traits than increasing (i.e. ‘winners’) or stable species did not receive support, as we did not find that species with certain traits would have become more common or rare in the study period.

Changes in the functional composition-environment relationships

Lake depth was the only variable explaining the spatial variation in functional community composition in both 1940s and 2017. It has been found earlier that the relative strength of competitive interactions in a specific functional niche could vary with water depth in macrophyte communities (Fu et al., 2014). In addition, colonisation depth and minimum light requirements vary with macrophyte growth forms (Middelboe & Markager, 1997), which was one of our functional traits. In addition to depth, the functional composition was associated with elevation and the amounts of ditches in the 1940s. In our study area, elevation represents lake landscape position and is temporally a very stable variable. The lake landscape position reflects both the connectivity and lake characteristics (Kratz et al., 1997; Riera et al., 2000), and in our study area, many environmental characteristics are related to this gradient. In other glaciated areas, macrophyte taxonomic community composition has been found to differ along the lake landscape position gradient (Alexander et al., 2008). Thus, it was not a surprise that this gradient was also important to functional composition of vascular macrophyte communities in our study area.

Ecke (2009) found that drainage ditching rather than land use itself could affect water quality and occurrences of macrophytes. However, it was interesting to note that the amounts of ditches were an important variable explicitly in the 1940s, as the amounts of ditches have increased towards 2017 (Fig. 3). It is thus possible that ditching has already had effects on macrophyte communities early on, and further ditching did not have additional impacts. Ditching can affect water chemistry and cause eutrophication in agricultural and urban areas, while in peatland areas it can influence water transparency (Ecke, 2009). In addition, ditches may facilitate the dispersal of aquatic macrophytes by acting as dispersal corridors (Soomers et al., 2010, 2013).

In contrast, in 2017, lake size and pH explained the spatial variation in functional community composition. Larger lakes sustain different functional composition than smaller lakes in our study. Thus, it can be assumed that larger lakes harbour contrasting habitats for biologically and ecologically different macrophyte species, thereby showing functional composition differing from that in small forest lakes. Similar studies on macrophyte communities do not exist to our knowledge, but findings on lake macroinvertebrates have shown that functional composition changes from small to large lakes along with concomitant changes in habitat structural features (Heino, 2008). The importance of pH is different from that of lake size. In acid and neutral waters, increasing pH might first enable the occurrence of different kinds of species causing increases in functional diversity, as high pH is related to increasing levels of bicarbonate and the availability of carbon for macrophytes (Vestergaard & Sand-Jensen, 2000; Iversen et al. 2019). Nonetheless, there is no clear spatial pattern in changes of pH from the 1940s to 2017 (Fig. 3), and based on the Wilcoxon signed rank test, there is no significant change in pH between the decades (Table 2). In addition, it should be taken into account that the pH measurements were done in different seasons, thus limiting the strict comparability between decades.

Patterns and drivers of functional richness

Zhang et al. (2018a) detected a decrease of functional richness of macrophytes in the Yangtze River floodplain, where changes in the environment have been drastic for a long time. We found a similar pattern in functional richness, even though changes in the environment have been only modest in our study area. We also found that the changes in functional richness through time were only partly linked to the concomitant changes in the environment, as stable characteristics of the environment—lake area, elevation and north coordinates—were more prominent than non-stable variables. As elevation in our study area presents the position of the lake in the landscape, it is not a surprise that it was collinear with other explanatory variables. Changes in functional richness were higher in the larger lakes, even though the pattern was not as clear as with elevation. One reason for this could be that, in our study area, the changes in human impact have also been more drastic in the larger lakes.

In our study area, the importance of north coordinate is most likely related to the glaciofluvial deposits and esker formation in the northern part of the study area and the historical development of the lakes nearby. These lakes nearby (e.g. lakes Kaukajärvi (10), Kirkkojärvi (11) and Tohloppi (25) in Fig. 2) are naturally mesotrophic clear-water and spring-fed to some extent. From the non-stable variables, only the changes in agricultural area had impact on changes in functional richness through time. Agriculture is linked to nutrient concentrations, especially total phosphorus (Ecke, 2009) and nitrogen-phosphorus ratios (Arbuckle & Downing, 2001). In our study area, the amount of agricultural area near the lakes has decreased through time, but the relationship with the change in functional richness and the change in agricultural area is not as clear as with elevation.

Traits of ‘losers’ and ‘winners’ through time

At the community level, anthropogenic impacts have been shown to increase the relative abundance of tolerant species with simultaneous decline or loss of specialist species (McKinney & Lockwood, 1999; Olden & Rooney, 2006). This would mean that species with certain sets of traits would become more common (‘winners’) or rare (‘losers’). However, at the species level, we did not find that this pattern had occurred in the study period. One reason for this could be a high intraspecific trait variation with aquatic macrophytes species, as many of them have a high phenotypic plasticity both morphologically and ecologically. A good example for this plasticity is Sagittaria sagittifolia L. (Lacoul & Freedman, 2006), which occurs in our study lakes. Due to plasticity in morphology and physiology, some species can have a different growth form in different environmental conditions (Wetzel, 2001). For example, Eleocharis acicularis (L.) Roem. & Schult. and Ranunculus reptans L. can occur as either helophytes or isoetids. Moreover, Persicaria amphibia (L.) Delarbre can grow either on land or in water (Hutchinson, 1975).

For macrophytes, Fu et al. (2014) found that intraspecific trait variability is an important factor influencing community dynamics and, in addition, intraspecific trait variation has been found to contribute to functional beta diversity patterns (Spasojevic et al., 2016). Therefore, it could be possible that macrophyte species can adapt or adjust to modest changes in the environment. In Germany, aquatic macrophyte species loss in streams has been associated with changes in species traits related to mechanical stress tolerance (Steffen et al., 2013), but the changes in the environment have been more drastic in Central Europe than in our study area in Northern Europe. On the other hand, if one considers the trophic status indicator features of these macrophyte species (Toivonen & Huttunen, 1995), most of the increased macrophytes prefer eutrophic and mesotrophic sites, indicating increasing eutrophication in time in lakes at lower elevations.

Why was functional homogenisation not detected in our study area?

Put together, there are probably four main reasons why we did not find signs of functional homogenisation in our study area during the 70-year study period. First, in our study area, where the changes in the environment have been quite modest, the stable factors in the landscape have a stronger effect than human induced changes on the vascular macrophyte functional composition and the changes in functional richness, thus masking functional homogenisation. Especially, the strong lake landscape position gradient is prominent in our study area. Elo et al. (2018) also did not find signs of biotic homogenisation in relatively oligotrophic lakes in Eastern Finland, where human impact is similarly quite low compared to many other areas worldwide. Second, changes in the environment in relation to human actions have been divergent through time. Urbanisation has intensified in our study area, while the agricultural land area near the lakes has decreased up to the 21st century (Table 2; Fig. 3). Thus, it is possible that these processes are masking the impacts of each other. Generally, in Finland, the use of fertilizers has also increased, while extensive areas of pastures on lakeshores has decreased and the municipal sewer systems have improved. Third, due to high phenotypic plasticity (Lacoul & Freedman, 2006), macrophyte species can probably adapt to these modest changes in the environment. Fourth, as primary producers, macrophytes can have greater landscape-level resilience than consumers (Johnson & Angeler, 2014). Johnson and Angeler (2014) found that macrophyte assemblages did not become homogenised with increased disturbance, although fish and macroinvertebrates did so. It could be argued that there is also a fifth reason. Jarzyna and Jetz (2018) found that temporal functional diversity patterns are scale dependent. Nevertheless, in lake environments, functional diversity patterns should be seen particularly at landscape level (e.g. Heino, 2008), which our study represents. Thus, we do not believe that the issue of spatial scale is a concern in our study. Furthermore, we consider that the time period covered by our study (~ 70 years) is sufficient to observe significant trends in aquatic macrophyte functional changes in case they would exist.

Conclusion

This study increases our knowledge on functional stability of vascular macrophyte communities and driving factors at the landscape level. Even though we did not find signs of functional homogenisation or differentiation, the changes in the environment have affected functional community composition and changes in functional richness to some extent. Thus, the effects of human activity on functional changes of macrophyte communities have probably become more prominent through time. As aquatic macrophytes have a crucial functional and structural role in lake environments, preserving diversity of functional features in aquatic plant communities should be taken into account in land use planning and conservation actions. It is also important to note that these results are based on four specific traits and other patterns could arise with other traits, even though the chosen traits are quite essential, robust and widely used in macrophyte studies. The selection of traits was still to some degree limited because there is a lack of trait information in many cases with more northern and oligotrophic species. Future studies on functional homogenisation should not concentrate only on areas where changes in the environment have been severe, but also on areas where human impact has not been that prominent to help us better understand functional homogenization process in its early stages.