No biotic homogenisation across decades but consistent effects of landscape position and pH on macrophyte communities in boreal lakes

294 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Thierry Oberdorff Editor-in-Chief: Miguel Araújo Accepted 25 September 2019 43: 294–305, 2020 doi: 10.1111/ecog.04757 doi: 10.1111/ecog.04757 43 294–305


Introduction
Biological communities are temporally dynamic, as they gain and lose species over time when the environments around them change (Bengtsson et al. 1997).
However, human actions have increased these changes by intensifying land use, introducing new species to ecosystems and altering the climate especially during the last century (Vitousek 1994, Chapin et al. 2000. One major consequence of human impact on biodiversity through time is that ecosystems are losing their biological uniqueness in a process called biotic homogenisation (McKinney andLockwood 1999, Olden andRooney 2006). McKinney and Lockwood (1999) first defined biotic homogenisation as 'the replacement of local biotas with non-indigenous species'. Subsequently, Olden and Rooney (2006) emphasised the multidimensional and multifaceted nature of this process occurring between two or more locations over a specified time interval, due to which ecosystems lose their biological uniqueness in space. Although biotic homogenisation has gained much attention in recent decades (Castaño-Sánchez et al. 2018, Richardson et al. 2018, it is still unclear how this process is acting at different levels of biodiversity through time in different ecosystems. This is especially the case in the lake-rich boreal region, which covers large areas across the Northern Hemisphere and is only moderately impacted by human activities compared to many other regions and ecosystems. One approach to tackle the issue of biotic homogenisation is to study changes in beta diversity patterns (Richardson et al. 2018). Beta diversity refers to the compositional dissimilarity in species assemblages across space or time (Whittaker 1972, Anderson et al. 2011. Many studies have focused on temporal beta diversity (i.e. within sites; Cook et al. 2018), while changes in spatial beta diversity (i.e. across sites) through time have gained less attention. The temporal patterns in spatial beta diversity can reveal how the similarity between sites changes over time (McGill et al. 2015). Therefore, beta diversity patterns can reveal the ongoing signs of degradation of ecosystems and the spatial aspects of biodiversity loss, biotic homogenisation (Rahel 2002, Olden andRooney 2006) or differentiation (Olden and Poff 2003). There is also growing evidence that beta diversity patterns generally could reveal more about changes in biodiversity compared with alpha diversity through time (Dornelas et al. 2014, Winegardner et al. 2017, Richardson et al. 2018. In addition, as contemporary and future threats as well as solutions in the conservation of natural environments occur at several spatial scales, beta diversity patterns can provide a useful tool in understanding these multidimensional issues (Socolar et al. 2016), especially in the world facing the threats of climate change and new species also spreading more rapidly to northern environments. McGill et al. (2015) predicted that spatial beta diversity (i.e. across sites) is showing a decreasing trend in the Anthropocene due to increasing human impact. However, they also noted that studies have shown that temporal patterns in spatial beta diversity have been highly context dependent (McGill et al. 2015). As these patterns are unclear, this issue should be studied more with different organism groups and in various environments.
Beta diversity patterns can be examined through different biodiversity facets, and recent studies have shown that using the taxon, phylogenetic and functional facets of beta diversity together can give better insights into spatial biodiversity patterns than the taxonomic approach alone (Heino and Tolonen 2017). In addition, biotic homogenisation is considered to cover the loss not only of taxonomic distinctiveness over time, but also of functional (Olden et al. 2004) and phylogenetic distinctiveness (Winter et al. 2009). This emphasises the fact that these different dimensions of biotic homogenisation should be explored simultaneously, even though most studies concerning biotic homogenisation have focused only on taxonomic homogenisation (Olden et al. 2018). Taxonomic distinctiveness refers to loss or replacement of native species, i.e. changes in beta diversity (Olden et al. 2004), while functional homogenisation refers to loss of specialised species or entire functional groups (Olden et al. 2004, Clavel et al. 2011. In phylogenetic homogenisation, endemic or rare species are lost, resulting in a decrease of among-species genetic differentiation (Winter et al. 2009). Thus, the use of different biodiversity facets (i.e. taxonomic, phylogenetic and functional) can reveal different aspects of biotic homogenisation and different mechanisms associated with these aspects of biodiversity.
It is estimated that changes in freshwater ecosystems globally, and especially in lake environments, will be due to biotic exchange, land use and climate change (Sala et al. 2000). However, many features of the physical environment of freshwaters can remain stagnant through time. For example, many characteristics of lake environments follow a pattern with the lake landscape position (Kratz et al. 1997, Riera et al. 2000, which rarely changes within short time periods. Position in the landscape can reflect both the hydrologic connectivity and physical features of a lake and the landscape, which are strongly related to lake water chemistry (Johnson et al. 1997, Martin andSoranno 2006). Thus, the lakes in one lake district can be very different regarding the abiotic characteristics and the biota. For example, the lakes lower in the landscape are typically larger, are more connected to the drainage system, and tend to have higher pH and ionic concentrations than the headwater lakes Muotka 2006, Martin andSoranno 2006). Therefore, it is possible that lake landscape position can overdrive, for example, changes in land use through time when beta diversity is considered.
We examined vascular aquatic macrophyte communities and spatial patterns in beta diversity in relation to human impacts from 1940s to 2010s. We had two main questions: Are spatial beta diversity patterns (i.e. across sites) different between decades, i.e. has there been biotic homogenisation? What are the main drivers for spatial beta diversity and are those drivers the same over the decades? We also aimed to discover if these patterns differ between species-, phylogenyand trait-based beta diversity (i.e. different beta diversity facets). We used presence-absence data of aquatic macrophytes from five different decades from small boreal lakes and used generalised dissimilarity modelling to study separately the relationship between different biodiversity facets and environmental gradients among five different time periods. We hypothesised that: 1) the patterns in beta diversity are different in relation to human impacts in different decades (Petsch 2016). As it has been predicted that human impact generally would cause a temporally decreasing trend in spatial beta diversity (McGill et al. 2015), we also hypothesised that 2a) the aquatic macrophyte communities have become more similar among lakes (biotic homogenisation), i.e. spatial beta diversity decreases due to human impact over time. As human impact has not been that severe in our study area, we alternatively hypothesised that 2b) the aquatic macrophyte communities have not become more similar among lakes, i.e. spatial beta diversity has not decreased due to human impact over time. 3) Human activities (i.e. land use changes) have a strong impact on beta diversity (Gámez-Virués et al. 2015). 4) The patterns in beta diversity are different for different facets of beta diversity (Heino and Tolonen 2017). Finally, 5) we also hypothesised that the strong landscape position gradient might overcome the effects of human land use change (Alexander et al. 2008) in driving beta diversity patterns. Studying these issues is important, as it has been found that spatial homogenisation of plant diversity reduces ecosystem multifunctionality (Hautier et al. 2018). In addition, freshwater ecosystems are globally considered biodiversity hotspots and provide valuable ecosystem goods and services to humans, but are intensively threatened by global change (Dudgeon et al. 2006, Vörösmarty et al. 2010, Vilmi et al. 2016, Reid et al. 2018).

Study area
The lakes studied were located in the Kokemäenjoki drainage basin in southern Finland near the city of Tampere, in the area between the two large lakes Roine and Pyhäjärvi (Fig. 1). The study area belongs to the southern boreal climate zone (Ahti et al. 1968), and it had a mean annual temperature of 4.4°C and a mean annual precipitation of 598 mm during the period 1981-2010 (Pirinen et al. 2012). The winter ice cover period lasted approximately 150-170 d during the period 1961-2000 (Korhonen 2005), and the length of the thermal growing season (> 5°C) was approximately 175-185 d during the period 1981-2010 (Finnish Meteorological Institute 2019). The underlying bedrock consists mainly of gneiss and diorites, soil mainly of sand moraine and rocks, and clay at lower elevations (Geological Survey of Finland 2018). The elevational gradient between the lakes studied is low (77-131 m). However, due to glacial and postglacial processes, the fine-grained sediment with nutrients was washed along the elevational gradient in the landscape (Seppälä 2005). Thus, changes in environmental conditions occur even along a relatively modest elevation gradient in the study area. Many of these lakes are situated in small chains of lakes and streams and have brown humid water, both features typical to Finnish lakes. More eutrophic lakes are located at lower elevations and are mainly surrounded by agricultural land or settlements. Smaller and more oligotrophic lakes at the higher elevations in the landscape are less affected by human activity and are mainly influenced by peatland drainage and summer cottages (Toivonen and Huttunen 1995). Key characteristics of these lakes are given in Table 1.

Species data
Aquatic macrophytes were sampled from 27 lakes using similar methods during five different decades. The first macrophyte surveys were conducted in 1947-1950by U. Perttula (unpubl.) and reinvestigated in 1975-1978(Toivonen and Huttunen 1995), 1991-1993-2008 For clarity, we refer to these surveys as 1940s, 1970s, 1990s, 2000s and 2010s. An aquascope and two different types of rakes were used in surveying aquatic macrophytes from the whole lake area. The macrophyte identification was done at the lowest possible taxonomic level. In this study, we concentrated on presence-absence data of species classified traditionally as aquatic vascular plants in Finland according to Linkola (1933) and, in addition, seven species growing in the water from the sedge genus Carex. In total, 68 vascular taxa were included (Supplementary material Appendix 1). We did not include hybrids and taxa identified to genus level in the analysis. We combined some species to species complexes due to identification differences between the five decades (Supplementary material Appendix 1). Species richness values and descriptive statistic in different decades can be found in Table 1.
In the absence of true phylogeny covering all the macrophytes in our data, we used taxonomic distances based on the Linnaean hierarchy as a proxy for phylogenetic relationship of macrophyte species. This approach has been used in previous studies dealing with phylogenetic diversity (Ruhí et al. 2013, Heino andTolonen 2017), but can only be considered as a coarse proxy of true phylogeny. We used equal branch lengths and five taxonomic levels above species level: genus, family, order, class and subdivision. Taxonomic information for macrophytes was collected from the open online source Catalogue of Life (Roskov et al. 2018).
To obtain functional dissimilarity of aquatic macrophytes, we used four biological traits: growth form, normal method of propagation, perennation and potential size. These are important biological traits of aquatic macrophytes (Willby et al. 2000, Göthe et al. 2017), affecting where the species can live, how they reproduce and what kind of life cycle they have. Growth form division consists of the following classes: ceratophyllid, elodeid, helophyte (incl. tall Carex), isoetid, lemnid and nymphaeid (Toivonen and Huttunen 1995). Normal method of propagation was collected from BiolFlor (Klotz et al. 2002). This trait consists of five classes: 1) by seed/spore, 2) mostly by seed/spore but also vegetatively, 3) by seed/spore and vegetatively, 4) mostly vegetatively and also by seed/spore and 5) vegetatively. Perennation information was mainly collected from Willby et al.'s (2000) attributebased data. 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 weighted 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 sources for all species. In such cases, data was complemented by information from other literature sources and databases. Perennation consisted of three classes: 1) annual, 2) biennial/short lived perennial and 3) perennial. Potential size information (cm) is a continuous trait (from Hämet- Ahti et al. 1998, Mossberg andStenberg 2012), representing the maximum potential length of an individual omitting the root or rhizome length (Bornette et al. 1994, Dolédec andStatzner 1994).

Environmental variables
We had only a limited amount of environmental information available from the 1940s, and thus we concentrated on environmental variables that are widely identified to be key variables for aquatic macrophytes: lake area (ha), maximum depth (m), pH and Secchi depth (m) (Lacoul and Freedman 2006, Table 1). All these environmental variables represent a larger complex of ecologically important factors. To study the effect of lake landscape position, we determined four different variables: elevation of the lake (m), watercourse distance to the main lake (m), lake order and lake network number. All these variables characterise the position of a lake along the watershed upland-lowland gradient and therefore several physical, hydrological and ecological characteristics of the lake (Kratz et al. 1997, Quinlan et al. 2003. All these four landscape position variables were correlated with each other based on Spearman's correlation test. Therefore, we did the modelling by using all these variables separately in our models. As the general patterns with each lake landscape variable were similar, we decided to focus only on elevation, as models that included it best explained variation in the compositional dissimilarity among lakes. In addition, we used land use variables from 200 m buffer zones (Pedersen et al. 2006) derived from the base maps for each decade (National land survey in Finland 2017, 2018) as a proxy for human impact. Several studies have shown that land use on relatively narrow buffer zones adjacent to lake shoreline has the strongest impact on lake plants (Pedersen et al. 2006, Akasaka et al. 2010. We calculated three proportional land use variables: agriculture area (i.e. field area), built area and ditches. We decided to focus on these land use types as they are key land use types that have changed most over the past decades in this study area. In addition, it has been shown that they influence water chemistry and other physical characteristics of lakes. The agricultural land within a watershed is strongly related, for example, to total phosphorus (Taranu and Gregory-Eaves 2008). Built area represents the human settlements and the level of urbanisation, which are often related to non-native species distribution (McKinney 2006). Ditches in the lake catchments have been shown to have an effect on lake macrophyte communities and water chemistry (Ecke 2009). In addition, we used geographical coordinates of lake centres. We also considered climatic variables; however, reliable climate data from our small study area were not available for the whole study period.

Statistical methods
A schematic diagram showing the methodology used can be found in Supplementary material Appendix 2. First, we calculated the amount of total beta diversity (reflecting both species replacement and loss/gain) among all pairwise comparisons of lakes (Legendre 2014) for taxon, phylogenetic and functional data and for different time periods separately. We calculated both the average and the variance of total beta diversity based on Jaccard's dissimilarity index (Legendre and Legendre 2012) using the function beta.multi in the 'BAT' package (Cardoso et al. 2015(Cardoso et al. , 2017. We decomposed the total beta diversity into replacement and richness difference components following the partitioning framework developed by Podani and Schmera (2011) and Carvalho et al. (2012). Because we were interested in overall difference in species richness, we applied the richness difference component (Podani and Schmera 2011) instead of the nestedness component (Baselga 2010) of beta diversity (Legendre 2014). However, we decided to focus only on total beta diversity, as the patterns in the replacement and richness difference components did not change clearly between the five decades. Second, we used generalised dissimilarity modelling (GDM) to analyse spatial patterns in different facets of total beta diversity in relation to environmental and geographical gradients (Ferrier et al. 2002(Ferrier et al. , 2007. GDM is a nonlinear extension of matrix regression for analysing and predicting patterns of compositional dissimilarity in relation to environmental gradients. It models compositional dissimilarity between all possible pairs of locations as a function of environmental differences between these locations. It takes into account the nonlinearity both in the relationship between ecological separation and observed compositional dissimilarity (Gauch 1973, Faith 1992 and in the rate of compositional turnover along environmental gradients (Ferrier et al. 2002(Ferrier et al. , 2007. To run GDM, we used three different pairwise dissimilarity matrices based on total beta diversity: dissimilarity matrix based on the aquatic macrophyte species presence-absence data (hereafter taxon dissimilarity matrix), phylogenetic dissimilarity matrix and functional dissimilarity matrix. We formed the taxon dissimilarity matrix based on Jaccard's dissimilarity index (Legendre and Legendre 2012) using the function beta in the 'BAT' package (Cardoso et al. 2015(Cardoso et al. , 2017. To obtain the phylogenetic dissimilarity matrix between sites, we first calculated taxonomic distances between aquatic macrophyte species using the function taxa2dist in the 'vegan' package (Oksanen et al. 2018). Then, we utilized hierarchical cluster analysis to produce a taxonomic tree for these species using the function hclust in the 'stats' package (R Core Team). Then we used this taxonomic tree along with the aquatic macrophyte absence-presence data and formed a phylogenetic dissimilarity matrix between sites based on Jaccard's dissimilarity index (Legendre and Legendre 2012) using the function beta in the 'BAT' package (Cardoso et al. 2015(Cardoso et al. , 2017. To obtain the functional dissimilarity matrix between sites, we first generated the species-by-species distance matrix based on the trait data using the function gowdis in the 'FD' package (Laliberté andLegendre 2010, Laliberté et al. 2014). This function is based on Gower (1971) (dis)similarity coefficient. Then, similarly to phylogenetic dissimilarity, we utilized hierarchical cluster analysis to produce a trait tree and subsequently generated the functional dissimilarity matrix based on the trait tree and aquatic macrophyte presence-absence data in the same way as the taxon and phylogenetic dissimilarity matrices.
Using the gdm package (Manion et al. 2018), we converted each dissimilarity matrix and the environmental data to site-pair format with the function formatsitepair. Then, we fitted the GDM model to tabular site-pair data using the function gdm and estimated the variable importance in the GDM model using the function gdm.varImp. The variable importance was quantified as the percent change in deviance explained by the full model and the deviance explained by a model fit with that variable permuted. We tested the full set of variables because we were particularly interested in the relative effects of these variables. This also facilitates comparing the impacts of different predictor variables with different beta diversity facets across all decades. In addition, we estimated the significance of each variable using the bootstrapped p-value using the function gdm.varImp. All statistical methods were conducted in R base ver. 3.4.3 (R Core Team), and the modelling was independently conducted for different time periods and different facets of beta diversity (i.e. taxon, phylogeny and functional).

Results
Gamma diversity has remained relatively stable in the study area between the decades. In the 1940s, there were 61 species, and the number of species increased by one per time point until the 2000s and in the 2010s there were 63 species. Alpha diversity has slightly increased from the 1940s until the 1990s, and has decreased after the 1990s (Supplementary material Appendix 5).

Multiple site beta diversity and components
There were no clear changes between the five decades in beta diversity. The amount of taxon beta diversity did not change much between the decades, as the average value ranged between 0.57 and 0.63 (Fig. 2, Supplementary material Appendix 3). The taxon beta diversity was driven by both the richness difference and the replacement. However, the richness difference component (65.7-68.5%) was more dominant than the replacement component (31.5-34.3%) in each decade (Supplementary material Appendix 3). The amount of phylogenetic beta diversity, on average (0.45-0.50), has not changed much between the five decades either (Fig. 2, Supplementary material Appendix 3), and it was also driven more by the richness difference (69.9-72.9%) than the replacement (27.1-30.1%) component (Supplementary material Appendix 3). Similar patterns occurred for functional beta diversity as the average values vary between 0.38 and 0.44 (Fig. 2, Supplementary material Appendix 3), and it was similarly more contributed by the richness difference (75.4-79.2%) than the replacement (20.8-24.6%) in each decade (Supplementary material Appendix 3). As the patterns in the replacement and richness difference components did not change clearly between the five decades and between the different facets of beta diversity, we concentrated only on total beta diversity in further analyses.

Predictors of different beta diversity facets
In the GDM models on taxon total beta diversity, the deviance explained by the environmental variables and geographical distances did not vary much between the different time periods. However, the explained deviance was clearly lower in the 1940s (56.4%) than in the later decades, when it was relatively high (76.4, 73.5, 79.5 and 69.5%, respectively) ( Table 2). In the 1940s, the largest amount of compositional change was observed along the gradient of elevation, followed by pH. In each following decade, pH was clearly more important than elevation, elevation still being the second most important (Fig. 3). Of the land use variables, only ditches in the 1940s and built area in the 2010s had a clear impact on taxon dissimilarity, but otherwise land use variables were not important. Geographical distance between lakes was a weak predictor of the compositional dissimilarity in each decade ( Explained deviance for phylogenetic total beta diversity models was only slightly lower than for the taxon-based GDM models in each decade. The explained deviance was lowest in the 1940s (49.3%) and highest in the 2000s (75.5%) ( Table 2). Elevation and pH were the two most important variables in each decade. As with the taxon-based models, elevation was more important than pH only in the 1940s (Fig. 3). Ditches were an important variable for phylogenetic dissimilarity in the 1940s. For phylogenetic dissimilarity, geographical distance between lakes was a poor predictor (Table 2, Supplementary material Appendix 4).
Explained deviance for functional total beta diversitybased models was lower in every decade than for the taxon and phylogenetic models. The explained deviance was lowest in the 2010s (42.0%) and second lowest in the 2000s (45.7%) and 1940s (47.5%) ( Table 2). pH and elevation followed the same patterns as the taxon and phylogenetic dissimilarity, except in the 2010s, when elevation was not highly important (Fig. 2). In the 2000s and 2010s, depth appeared as an important variable. Of the land use variable, ditches in the 1940s and built area in the 2010s had impacts on functional compositional dissimilarity, as was the case also with the taxon and phylogeny data. For functional dissimilarity, geographical distance between lakes was again a poor predictor. As a whole, the predictor variables' impacts on functional dissimilarity diverged slightly from the other two facets of beta diversity (

Discussion
It has been predicted that spatial beta diversity (i.e. across sites) is showing a decreasing trend in the Anthropocene due to human impact (McGill et al. 2015). However, we did not find signs of such decreasing trend in our study area during the past 70 yr. Instead, we found that vascular aquatic macrophyte communities showed only slightly different patterns in spatial beta diversity across decades, thus contradicting H1. Many studies have shown that human actions are causing biotic homogenisation in different environments in general (Richardson et al. 2018) and in freshwaters in particular (Castaño-Sánchez et al. 2018, Zhang et al. 2018b), and Salgado et al. (2018 found this trend also for aquatic macrophytes in shallow lakes. Nevertheless, we did not find signs of biotic homogenisation (Olden and Rooney 2006) or biotic differentiation (Olden and Poff 2003). This observation suggests that changes in land use and other human impacts have not homogenised these communities at the landscape level (contradicting H2a, supporting H2b). Moreover, we did not find changes in species relatedness or functional similarity of communities across decades, suggesting that neither phylogenetic (Winter et al. 2009) nor functional homogenisation (Clavel et al. 2011) are occurring in our study area.
In addition to beta diversity, there has not been clear changes in gamma diversity in the study area between the decades. However, the median species richness (alpha diversity) has increased until the 1990s and has decreased after Figure 2. The average of the total beta diversity for the taxon, phylogenetic and functional dissimilarity between the sites in each decade. 301 that decade. As the species richness difference was more dominant than the replacement component in each decade and for each beta diversity facet, it suggests that there is a strong species richness gradient that persists across decades, despite the changes in local species richness and the environment. The temporally persistent gradient is associated with lake landscape position, with lakes high in the landscape having low species richness and lakes low in the landscape have generally high species richness. This pattern is likely to result from the combined effects of local environmental conditions and dispersal on macrophyte species richness.
Studies that have found clear signs of biotic homogenisation for aquatic macrophytes are either considering other biodiversity measurements or beta diversity measures, or are based on paleoecological methods (Salgado et al. 2018). There are only a few studies of temporal changes in spatial beta diversity (Winegardner et al. 2017, Larsen et al. 2018, Wengrat et al. 2018) and even fewer considering aquatic macrophytes. McGill et al. (2015) proposed that the patterns of spatial beta diversity were unclear and that there were strong influences of context dependency among studies. For example, Winegardner et al. (2017) did not find changes in spatial beta diversity of lake diatoms between 1850 and 2007, even though they had a longer time period examined compared with that in our study. Similarly, Larsen et al. (2018)   did not find signs of biotic homogenisation when they studied stream macroinvertebrates in 10 streams during 30 yr. However, Wengrat et al. (2018) found a decreasing trend in spatial beta diversity of diatom assemblages, but only when they studied eutrophic reservoirs instead of the whole set of reservoirs over the past 60-100 yr.
Compared to the other areas in the world where biotic homogenisation has been observed (Zhang et al. 2018b, Finderup Nielsen et al. 2019, the environment in our study area has faced only modest changes from the 1940s to the 2010s. Admittedly, we found that land use was not particularly important in explaining the compositional dissimilarity in our study area, thereby contradicting H3. During the last 70 yr, the agricultural land area has decreased and simultaneously urbanisation has intensified in our study area (Table 1). Thus, it is possible that the opposite effects of these factors could balance each other in affecting beta diversity. Both land use types are known to increase pollutants and nutrients in freshwaters, further affecting biota (Paul andMeyer 2008, Taranu andGregory-Eaves 2008). The built area, representing the human settlements and the level of urbanisation, has even showed a small increase in importance in explaining compositional dissimilarity after 1990s. Also, Elo et al. (2018) did not find biotic homogenisation in relatively oligotrophic lakes in eastern Finland, where human impact is also low compared to other areas globally. On the other hand, Zhang et al. (2018a) found taxonomic and functional differentiation of macrophyte assemblages instead of homogenisation, even though their study area in the Yangtze River floodplain has faced strong human impact for a long time.
It is also possible that relatively modest changes in land use do not result in altered beta diversity patterns based on presence-absence data. Using only binary coefficients instead of quantitative forms of the indices can also produce coarser results (Legendre 2014). Even though we did not find changes in spatial beta diversity patterns, it is possible that there have been changes in species abundances. For example, there could be a loss of functional beta diversity due to more specialized species becoming rare but not extinct in the lakes. However, when using 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 also hypothesised that the strong landscape position gradient might overcome the effects of these quite modest land use changes in driving beta diversity patterns (H5). We used elevation to represent the landscape position, and it was an important variable explaining the compositional dissimilarities of all facets, especially in the 1940s. Thus, lake position in the landscape, which reflects both the connectivity and the lake characteristics (Kratz et al. 1997, Riera et al. 2000, probably explains partly the patterns found in beta diversity, thus supporting H5. In our study area, the more eutrophic lakes are situated at lower elevations, and they are surrounded by human settlements or arable lands. In contrast, the more oligotrophic lakes are situated in the upper parts of the drainage basins, and they are surrounded by coniferous forest and peatlands with less human impact. Moreover, an earlier study based on partly the same macrophyte data in the 1970s showed that nutrient concentrations and macrophyte species richness increased with the eutrophy of lakes along one of the lake chains (Ilmavirta and Toivonen 1986). Generally, lakes higher in the landscape are smaller, less connected, have lower stream inputs Muotka 2006, Martin andSoranno 2006; Fig. 1), and have lower species richness than the lakes lower in the drainage basin (Kratz et al. 1997). It has been previously found that macrophyte community composition differs along the lake landscape position gradient (Alexander et al. 2008). Additionally, the dominance of the species richness difference gradient also supports the lake landscape position hypothesis (H5). However, it was also notable that the importance of elevation decreased through time.
For all facets of beta diversity, pH was the most important variable in each decade, except in the 1940s when elevation was more important. Several studies have shown that pH has a major influence on aquatic macrophyte communities (reviewed by Lacoul and Freedman 2006). Thus, it was not surprising that it was an important predictor variable. Nevertheless, pH had a lower influence on beta diversity in the 1940s and also in the 2000s. Landscape characteristics, especially bedrock and soil, have a strong impact on pH values of lakes (Quinlan et al. 2003). In our study area, soil is linked to the lake landscape position as the clay areas are more or less located at lower elevations and the sand moraine and rocks are situated in higher elevations (Geological Survey of Finland 2018). In addition to pH, other water quality metrics, for example nutrients, play a prominent role for aquatic macrophytes (Lacoul and Freedman 2006). Thus, these findings may be somewhat limited by the fact that we did not have nutrient information available (from the 1940s). However, as the land use and nutrient concentrations are closely related (Johnson et al. 1997, Taranu and Gregory-Eaves 2008, Ecke 2009), this should not be a critical issue in our study.
It was interesting that the patterns related to different facets of beta diversity diverged only slightly from each other, which was mostly incongruent with H4. Taxon and phylogenetic facets of beta diversity followed the same patterns overall. Most notable differences were found for functional beta diversity, even though those differences between functional and taxon or phylogenetic beta diversity were quite small as well. Also, Zhang et al. (2018a) found that functional diversity did not add any particular value compared to taxonomic diversity. However, studies focussing on other organisms have shown that even using these kinds of quite coarse proxies of phylogenetic and trait information can provide additional information to taxon-based information Tolonen 2017, Richardson et al. 2018). Naturally, using true phylogeny could provide a different picture from taxonomic distances based on the Linnaean hierarchy. Moreover, many aquatic macrophyte species have high phenotypic plasticity (e.g. Sagittaria sagittifolia, Lacoul and Freedman 2006) and, additionally, intraspecific trait variation has been found to contribute to functional beta diversity patterns (Spasojevic et al. 2016). However, the traits we used were quite robust related to these issues studied, as we used quite broad trait classes and size categorization, specifically the potential maximum size. In addition, such trait divisions have been used repeatedly with aquatic macrophytes to represent functional diversity (Zhang et al. 2018a). Nevertheless, it has been found that temporal taxonomic and functional diversity patterns are scale dependent (Jarzyna and Jetz 2018). Thus, it is possible that changes in functional beta diversity patterns do not appear at the spatial and temporal scales examined in our study.
By using the spatial beta diversity perspective, our study highlights the fact that even though biotic homogenisation is a pervasive problem globally (Olden et al. 2004), it is not an unambiguous process acting similarly at all spatial and temporal scales or in different environments and different organism groups. When comparing our results to other studies, it is clear that when studying beta diversity patterns and biotic homogenisation, the context dependency, the degree of human pressures, the scale of the study (both spatial and temporal) as well as the measurements of biotic homogenisation are important when results are interpreted. Our findings further suggest that lakes across the boreal region and areas that have faced glaciation and postglacial processes might be resistant against moderate levels of human pressure. However, is also important to emphasize that our present results do not imply that there could not be biotic homogenisation at an individual lake level or across a longer time period. There might still be changes in temporal beta diversity patterns in individual lakes (Winegardner et al. 2017), and this issue should be further explored in boreal lake environments. Further studies should also take into account the lake landscape position and the issue of connectivity to increase our understanding of the landscape-scale biodiversity patterns. This study is highly unique in lake environments as the temporal aspect is often neglected in aquatic biodiversity studies. More specifically, our findings help understand how aquatic macrophyte communities may respond to changes in the environment across decades. Although we focused on a single lake district, the patterns detected in macrophyte beta diversity within and across decades are likely to represent situations in the large boreal areas of Eurasia and North America.