Shifts in Lichen Species and Functional Diversity in a Primeval Forest Ecosystem as a Response to Environmental Changes

Research highlights: shifts in the composition and functional diversity of lichen biota reflect changes in the environment caused by climate warming and eutrophication. Background and objectives: studies on lichen functional diversity and refinement in the functional traits of lichen biota under the pressure of changing environmental factors are currently of great scientific interest. The obtained results are interpreted in relation to specific habitat properties and their modifications due to the potential effects of climate change and atmospheric pollution. The aim of the work was to investigate changes in lichen species composition and functional diversity, as well as to identify factors responsible for them at different forest ecosystem scales. Materials and Methods: we identified factors responsible for changes in lichen biota in a unique Białowieża Forest ecosystem by analyzing shifts in species optima and functional diversity at the forest community, tree phorophyte, and substrate levels. We examined individual lichen species’ responses and temporal shifts in the species composition for each historical and resampled dataset using a communityweighted means of functional lichen traits and Wirth ecological indicator values. Results: the most evident change took place at the level of individual species, which shifted their realized optima: 25 species demonstrated a shift to co-occur with lichens of higher nitrogen demands, 15 demonstrated higher light demands, 14 demonstrated higher temperature preferences, and six demonstrated lower moisture preferences. At the level of forest communities, biota shifted towards the higher proportion of nitrogen-demanding and the lower proportion of moisture-demanding species. At the level of phorophyte species, biota changed towards an increased proportion of lichens of higher temperature preferences. For the substrate level, no directional shifts in lichen species composition were found. Conclusions: climate change has influenced lichen biota in Białowieża Forest, but the main driver of lichen species composition was found to be eutrophication. We suppose that other overlapping factors may contribute to biota shifts, e.g., the extinction and expansion of phorophyte tree species.


Introduction
Lichens are extremely vulnerable to any modification in their environment, to which they react by changes, for example, in species composition and abundance [1]. Among the various factors shaping lichen diversity, land-use and habitat transformation [2][3][4], forest fragmentation [5], biological invasions [6,7], urbanization [8], and climate change [9] are considered the most important. Accelerating climate change [10] is the main driver that is shaping lichen biota in many geographical areas [1,9,11,12]. Numerous lichenological reports have described transformations of biota caused by increasing CO 2 emissions, rising temperature, and decreasing precipitation [1,11,12]. Additionally, the progressive eutrophication of habitats due to the deposition of nitrogen compounds plays an important role in environmental conditions in which the species occur now have shifted over time compared to environmental conditions in the past. This indicates transformations in the environment, and it also provides information on how species' realized niches and co-occurrence patterns, as well as the whole composition of lichen communities, may have changed [21][22][23].
An earlier paper by Łubek et al. [14] on Białowieża Forest revealed shifts in the species composition of epiphytic lichens towards an increasing contribution of species with higher nitrogen requirements, which was related to an increasing deposition of nitrogen pollutions and the associated process of habitat eutrophication [25][26][27]. Despite the significant increase in air temperature (climate warming) and decrease in precipitation (habitat drying) over the last 30-50 years [26][27][28], modifications in epiphytic lichen biota were not associated with these climatic parameters. The effects of climate change on lichen biota are often not detectable in the case of the simultaneous reaction of all species, especially if these shifts are at the initial stage and very gentle [14]. Therefore, there is a need to separately investigate the response of each species and changes in its frequency using an individual approach, as species with broad habitat preferences may benefit from the transformation of environment and colonize new habitats, while specialists may reduce their frequency or even disappear [29]. Therefore, it is necessary to study the environmental optima of the historical and present co-occurring species at the level of individual species, as has been done for plants [21][22][23], which may reflect the species shifts under the changing habitat properties.
In our work, we focused on changes in the species optima and total lichen functional diversity (including epiphytic, epixylic, and epigeic lichens) in different forest communities and on different phorophytes and substrates to identify developments in lichen biota species composition and important factors that may potentially alter lichen composition over time. The aim of our work was: (i) to investigate changes and identify factors causing them at the ecosystem scale by analyzing shifts in species' realized optima; (ii) to identify forest communities, phorophytes, and substrates with the greatest adjustment; (iii) to explore whether changes at the species level correspond to those at the community, phorophyte, and/or substrate levels.

Study Area
The study area is located in the north-eastern part of Poland in Białowieża Forest, which is the best preserved forest area on the European lowland [30][31][32][33]. Białowieża Forest is characterized by high structural and spatial forest complexity, as well as the continuity of biological processes, so it provides suitable habitat conditions for a huge diversity and species richness of plants, animals, and fungi. The best preserved forest communities of Białowieża Forest are protected in Białowieża National Park (BNP), excluded from direct human impact since 1921. In BNP, various tree species reach a considerable size, dying and dead trees are the source of a large amount of dead wood, and fallen trees and the exposed roots of windthrows shape specific microhabitats not found in managed forests. All of these components, as well as the high moisture content of decaying wood, create excellent habitats for lichens. The high degree of the forest's naturalness is reflected in the high richness of lichen species and the presence of many lichen relicts of the primeval forests [34][35][36]. The area of Białowieża Forest, due to a continuity of forest cover [37] and a lack of direct human activity, is excellent for lichenological studies [7,14,17]. Białowieża Forest is characterized by a transitional temperate climate-it is relatively cool, with continental influence, and it is defined as a forest of subcontinental climate of the temperate zone [38]. Climate change has already affected the Białowieża Forest. Published data clearly indicate a decrease in the number of days with the minimum air temperatures and an increase in the number of days with the maximum temperatures over the years of 1950-2007, which has been expressed in the increase of the mean annual air temperature by 1.1 • C (from 6.3 to 7.4 • C) [26][27][28]39]. Changes have also been observed in the amount of precipitation, which was lower in the period of 1986-2007 (606 mm) compared to    [26,27,40]. The duration of snow cover has also significantly reduced by about 30 days [28,40]. In the case of atmospheric pollutants, decreases in sulphur deposition from approximately 15 kg ha −1 year −1 (dry deposition in gaseous form) and approximately 12 kg ha −1 year −1 (wet deposition in precipitation) in 1986 to approximately 1 and 6 kg ha −1 year −1 in 2007, respectively, have been observed. In the case of nitrogen, the dry deposition has fluctuated around the level of 2-3 kg ha −1 year −1 , and the wet deposition increased from 7-9 kg ha −1 year −1 in 1986-1995 to 14 kg ha −1 year −1 in 2005 [26,27]. This may be of particular importance for the eutrophication of habitats in Białowieża Forest [41], which is additionally amplified by habitat drying [26,27].

Data Collection
The field studies were conducted in BNP (52 • 46 N 23 • 52 E) on a set of 144 permanent plots, 100 × 100 m each [42]. In the investigated area, six of the most common forest communities representative of the whole of BNP occurred: mixed deciduous fertile oaklinden-hornbeam forest (54 plots), floodplain streamside alder-ash forest (22), swamp alder carr (18), coniferous pine-oak mixed forest (27), coniferous moist oak-spruce forest (18) and coniferous mesic (spruce)-pine forest (5). Of the 144 plots, 64 were covered entirely by only a single community. In the remaining plots, 2-3 different communities occurred, but one forest community dominated and covered at least 70% of the plot. Some of the communities were found to be similar in spatial structure, i.e., they were composed of similar tree species, contained similar types and amounts of substrates, and were similar in terms of microclimate conditions [14,17,43,44]. These similarities influence the grouping of lichen species of specific functional traits [14,17,45]. Therefore, in order to obtain clear results, we combined the investigated forest communities and put them into three groups: mixed deciduous forests (dominated by a mixture of deciduous tree species; 54 plots), wet forests (comprising floodplain streamside alder-ash and swamp alder carr; dominated by black alder and European ash; 40 plots,) and coniferous forests (comprising a coniferous pine-oak mixed forest, a coniferous moist oak-spruce forest, and a coniferous mesic (spruce)-pine forest; dominated by Scots pine and Norway spruce; 50 plots [17]).
The original survey was conducted in 1987-1989 [45], and a resurvey took place in 2014-2015. In both surveys, the same research method was used: in each plot, all lichen species (epiphytic, epixylic, and epigeic) on all available substrates were recorded [14,42]. We analyzed lichen biota on the most frequent phorophyte species: European hornbeam Carpinus betulus L., small-leaved lime Tilia cordata Mill., European ash Fraxinus excelsior L., Norway maple Acer platanoides L., European white birch Betula pendula Roth, European oak Quercus robur L., European aspen Populus tremula L., European alder Alnus glutinosa Gaertn., common hazel Corylus avellana L., Norway spruce Picea abies (L.) H. Karst, and Scots pine Pinus sylvestris L. For epiphytic and epixylic lichens, diverse types of substrates were analyzed: the bark of trunks of living trees (from 40 to 200 cm in height), the foot of living trees (up to 40 cm in height), the branches of living trees (all branches available up to 200 cm in height), the stems and branches of shrubs, the bark of trunks and the branches of fallen trees, the wood of trunks and the branches of fallen trees, the wood of dead standing trees, small tree stumps, and the root systems of fallen spruces. In the case of substrates for epigeic lichens, only humus and conifer litter in coniferous forests were available. The frequency of species was measured as the number of plots on which the species were recorded.
For species identification, a hand lens with a light was used in the field. For species requiring anatomical and chemical (thin layer chromatography) analyses for their proper determination, samples were collected for further examination (e.g., Ochrolechia, Cetrelia, Cladonia, Lepraria, and Micarea species, as well as many sterile lichens). During the resurvey, we found dozens of new species that had not been previously reported from the study area. Many of them were not known to lichenologists before 1990, as these species were described during the last three decades (e.g., Biatora pontica Printzen & Tønsberg, Lepraria jackii Tønsberg, and Micarea viridileprosa Coppins & Van den Boom) or were not recognized by previous surveyors due to the lack of the techniques used widely at present in lichen determination, such as thin layer chromatography (e.g., the case of the Cetrelia monachorum (Zahlbr.) Culb. & C. Culb. and Cladonia chlorophaea (Flörke ex Sommerf.) Sprengel groups). These species were certainly present in the study area during the original survey, but due to the lack of complete information on their occurrence, as well as to reduce biases linked to the observer effect on the interpretation of results, they were removed from analyses or combined into sensu lato groups (e.g., Cetrelia cetrarioides s.l., Cladonia chlorophaea s.l., Lepraria incana s.l., and Micarea prasina s.l. [14]).
The nomenclature of species mostly follows the work of Fałtynowicz and Kossowska [46]. The collected material was deposited in the lichen herbaria of the Jan Kochanowski University in Kielce (KTC) and the University of Gdańsk (UGDA).
In the analyses, we also used information about the ecological requirements of lichen species obtained from Wirth (  (Table S1; compare [14,17]. These indicators use a nine-step scale, with 1 being the lowest and 9 the highest values of the respective environmental factor, and the values can be indirectly used as additional functional traits that provide information on the range of lichen species requirements for a particular environmental factor.

Species Optima along Environmental Gradients
Due to the lack of plot-specific data describing environmental conditions in the historical dataset, we examined the individual species' responses to the underlying shifts in the environment using an indirect approach-the so-called 'weighted average technique' (sensu Kapfer et al. [21] Kapfer and Popova [23]). Using the weighted average technique, we explored relative shifts in lichen species' realized optima along different environmental gradients, represented by light, temperature, continentality, moisture, and nitrogen WEIVs. Because species individually respond to changes in the environment, part of the information on potential drivers that may shape lichen composition along environmental gradients may be lost when considering site-specific scales [53]. Therefore, we independently conducted this analysis for forest communities, tree phorophytes, and substrates. To avoid obtaining misleading results provided by rare species, we only tested those lichens that occurred in at least five plots in both the initial and resurveyed datasets, reducing the number of analyzed species from 177 to 102.
Considering all lichen species in a plot (with the exception of the focal species for which the optimum was being computed) for each period of data collection (historical vs. resampled) and for each WEIV, we computed a weighted average sample score. To ensure an equal distribution of sample scores along particular WEIV gradients, we removed extreme values of samples distributed outside the common range of both historical and resampled datasets, thus obtaining two datasets with equal maximum and minimum values for each weighted average sample score. Next, we ordered all weighted average sample scores into three equal categories along each WEIV gradient. In each category, we equalized the numbers of plots between the samples representing two periods of sampling with the most samples, thus obtaining two pruned datasets with similar ranges and frequency distributions of samples along each WEIV gradient. This procedure allowed us to compute a weighted average species score for each species (hereinafter: 'species optimum'-SOV), with a restricted permutation test performed to examine whether obtained shifts in SOVs were higher than expected by chance (with a p-value level at 0.05). A positive shift in SOV may have indicated that a species recently co-occurred more often with species of a higher value of specific ecological indicator or functional trait than historically, while a negative shift in SOV may have expressed the opposite trend. Thus, if shifts in SOVs in respect to one particular gradient were directional (i.e., most of the species revealed a positive or negative shift in SOV), a change in the ecosystem along this environmental gradient was also directional, i.e., towards higher or lower values of the analyzed parameter [21]. The same weighted average technique was employed to discover shifts in species' optima along the biotic gradients of lichens' functional trait values (only numerical traits we included in this procedure), i.e., ascomata area, ascospore shape, and ascospore volume, thus allowing us to make more complex estimations of shifts in the species' realized niches.

Species and Community Characteristics
To explore temporal shifts in the role of niche partitioning and habitat filtering in shaping the species composition of lichen biota, for each historical and resampled plot representing each forest community, each tree phorophyte, and each substrate we calculated community-weighted means (CWMs) of functional lichen traits (only numerical traits were included in this procedure) and Wirth ecological indicator values (considered as additional information on ecological requirements of lichen species describing their realized niche). In addition, using the entire set of functional traits (containing both categorical and numerical traits), we calculated four parameters of functional diversity (FRic, FEve, FDiv, and FDis) using the FD;:dbFD function [54]. Functional richness (FRic), evenness (FEve) and divergence (FDiv) were computed according to the method of Villéger et al. [55], while functional dispersion was calculated according to the work of Laliberté and Legendre [56]. These parameters describe different aspects of the distribution of functional traits in the community trait hyperspace [56]. While FRic is a metric of the filled niche size [57], FEve describes the regularity of trait distribution in a niche space [58]. FDiv and FDis are metrics of distances between values of different functional traits belonging to different species in a community hypervolume [57]. FDiv, in comparison to FDis, provides more detailed information on the role of specialized species with extreme values of functional traits in a community [17,59]. High values of functional diversity parameters can indicate a high richness and divergence of life strategies realized by different lichen species, thus expressing the high importance of niche partitioning and the low role of habitat filtering (environmental stress) in shaping community composition [60]. Low values of functional diversity indices can be an expression of the high convergence of life strategies realized by different species under environmental stress conditions (high habitat filtering and low niche partitioning), indicating that one or a few functional types of lichens are prevalent in a community [61]. To assess the taxonomic diversity of lichen species in the study site, for each plot representing each forest community, as well as for each tree phorophyte and substrate, we calculated the species richness and Shannon index.

Species Composition at Community, Tree Phorophyte, and Substrate Levels
Using ordination techniques, we explored temporal shifts in lichen species composition in more detail by considering the forest community, tree phorophyte, and substrate levels separately. The preliminary detrended correspondence analysis (DCA; vegan::decorana() function; [62]) performed for each forest community revealed relatively short gradient lengths in each case (<2.5 SD units). Therefore, we decided to use PCA ordination (vegan::rda() function) as the most proper exploratory method in assessments of temporal shifts in lichen species composition along environmental gradients at the forest community level. Due to the relatively long gradient length (>2.5 SD units) obtained for both the tree phorophyte and substrate levels, we chose DCA ordination as the most appropriate method for exploring changes in the species composition of lichens over time regarding both levels of ecosystem organization. Due to the low ecological importance of results driven by the inadequate distribution of sites (e.g., tree phorophytes and substrates) and species in an ordination space, prior to DCAs, we transformed the data into a binary form (vegan::decostand() function). To assess the indirect relationships between the distribution of study sites (plots representing forest communities, tree phorophytes, and substrates) and sampling period, lichen species taxonomic (species richness and Shannon index), functional diversity (FRic, FDiv, FDis, and FEve), CWMs of WEIVs (WEIV.L, WEIV.T, WEIV.C, WEIV.M, and WEIV.N), and lichen functional traits (only numerical traits: ascomata area, ascospore shape, and ascospore volume), we applied the vegan::envfit() function. Thus, we performed a passive projection of the above-mentioned variables fitted as vectors into an ordination space. Using a permutation test with 999 iterations, we computed the determination R 2 coefficient and p-value for each variable.
To compare values of the Shannon index and parameters of functional diversity, as well as to check for potential changes in the CWMs of Wirth indicators and lichen functional traits (only numerical traits we included in this procedure) between the two time periods of sampling at the forest community level, we used linear mixed effect regression with Gaussian distribution (LMM; lme4::lmer() function [63]). To assess changes in plot species richness over time, we employed generalized linear mixed effect regression with Poisson distribution (GLMM; lme4::glmer() function). Thus, for each analyzed parameter regarding each forest community, we performed a separate LMM (or GLMM in the case of species richness), with the sampling period adopted as a fixed effect and plot identity adopted as a random effect, to include plot-specific factors that could potentially influence the interpretation of results. These plot-specific factors accounted for the effects of the forest community constituting the minority of coverage in plots dominated by one forest community reaching at least 70% of the cover and effects of edge zones, i.e., the neighborhood of one plot with a particular community by plots representing different forest types. To check how much of variance was explained by the sampling period alone (fixed effect), for each LMM (or GLMM), we calculated the marginal coefficient of determination (R 2 m ). To reveal how much of variance was explained by both sampling period and plot identity (random effect), we computed the conditional coefficient of determination (R 2 c ). Thus, the difference between R 2 c and R 2 m indicated how much of variance was explained by the random effect alone (MuMIn:;r.squaredGLMM() function [64]). Due to the lack of data on the tree phorophyte, substrate diversity, and distribution for each plot, we did not perform LMMs (or GLMMs) in order to discover temporal shifts in lichen species composition at the tree phorophyte and substrate levels. All analyses were performed in the R software [65].

Shifts in the Realized Species' Optima
Of the 102 lichen species analyzed for relative shift in SOV, 52 were found to co-occur with lichens of different WEIVs compared with the initial sampling period (Figure 1a). Of the 20 species that demonstrated a significant change in SOV for light WEIV, 15 species shifted towards higher values for this SOV. Out of the 16 lichens that revealed a significant change in SOV along the gradient of temperature WEIV, 14 species co-occurred more often with more warm-demanding species. Regarding the continentality gradient, of the 17 lichen species that significantly changed their SOV for this WEIV, 15 species were found to co-occur with lichens of lower continentality preferences. Of the nine species that significantly shifted their SOV for moisture, six lichens were found to co-occur with lichens of lower values for moisture WEIV. For all 25 lichen species that had a significant change in SOV for nitrogen, we observed an increase in co-occurring lichens with a higher WEIV.N.
Considering SOV for lichens' functional traits, we found that in 1987-1988 and 2014-2015, the species' optimum significantly changed for 38 lichens (Figure 1b). Of the 12 lichen species that revealed a significant shift in SOV for ascomata area, 11 species co-occurred with lichens with higher values of this trait. Regarding the ascospore shape gradient, of the 15 lichen species that changed significantly their SOV for this trait, 12 lichens were found to co-occur with lichens of lower values of this parameter. We also observed that 16 species significantly changed their SOV for ascospore volume, and 13 of them were recently found to co-occur more often with lichen species of higher values of this parameter.

Community-Specific Changes in Species Composition
The principal component analyses performed for each forest community revealed large differences in the lichen species composition over time. We found the main shift in species composition between the two time periods of lichen biota sampling along PCA axis 1 in the case of mixed deciduous and wet forest communities, while this shift took place along both PCA axes for plots representing coniferous forests.
In all forest communities, the resampled species composition of lichens was characterized by higher values of vectors representing species richness, the Shannon index and FRic, as well as lower values of FEve (Figure 2a-c; Table S2). In plots representing mixed deciduous ( Figure 2b; Table S2) and wet forest communities ( Figure 2c; Table S2), these two last variables demonstrated a strong negative correlation with each other, while we observed no correlation between FRic and FEve for coniferous forests (Figure 2a; Table S2). In the same forest types, higher values of the vector representing FDiv were found in the recent survey (Figure 2a; Table S2) Table S2). We also found that the species composition of some resampled plots representing wet forests type was characterized by higher values of the vector representing light WEIV (Figure 2c; Table S2). Considering vectors representing WEIVs for temperature and continentality, as well as the CMWs of ascomata area, ascospore shape, and ascospore volume, we only found weak relationships with the lichen species composition of the historical and resampled plots (Figure 2a-c; Table S2). S2). We also found that the species composition of some resampled plots representing wet forests type was characterized by higher values of the vector representing light WEIV (Figure 2c; Table S2). Considering vectors representing WEIVs for temperature and continentality, as well as the CMWs of ascomata area, ascospore shape, and ascospore volume, we only found weak relationships with the lichen species composition of the historical and resampled plots (Figure 2a-c; Table S2).  Table S1 in the Supplementary Materials.  Wirth's ecological indicator values) of light, temperature, continentality, moisture, and nitrogen (a) and biotic environmental gradients (calculated based on lichen functional traits) of ascomata area, ascospore volume, and ascospore shape (b). Numbers in bars represent the size of the change in species optimum. For full species names, see Table S1 in the Supplementary Materials. In all forest communities, in comparison to the initial sampling period, we recorded significantly higher species richness, Shannon index, FRic, and the CWM of nitrogen WEIV values, while we reported significantly lower values for FEve and the CWM of moisture WEIV (Figure 3; Tables S3-S5). We identified the wet forest type as the community in which we found the highest amounts of variance explained by the sampling period (fixed effect; ranging from 43% to 79%), and the lowest amounts of variance explained by plot identity (random factor; ranging from 2% to 8%) in the case of almost all parameters mentioned above (except for FEve and the CWM of moisture WEIV; Tables S3-S5). Considering FDis and FDiv, significant increases in their values over time were only recorded in plots representing coniferous communities, with 9% and 10% of variance explained by sampling period and 69% and 3% of variance explained by plot identity, respectively ( Figure 3; Table S3). The only community for which we recorded a significant increase in the values of light WEIV CWM was the wet forests type, with 13% of variance explained by sampling period and 23% of variance explained by plot identity (Figure 3; Table S5). Regarding the CWM of ascospore shape, we found significant increases of its values over time in the wet and coniferous forests. About 12% and 11% of variance of this parameter were explained by the fixed effect and 21% and 34% were explained by the random factor in wet and coniferous forests, respectively (Figure 3; Tables S3 and S5). Additionally, in resampled plots representing coniferous forest communities, we observed significantly higher values of ascospore volume CWM, with 4% of variance explained by fixed effect and 24% of variance explained by random factor (Figure 3; Table S3).

Shifts in Species Composition at Tree Phorophyte and Substrate Levels
The detrended correspondence analysis performed for tree phorophytes revealed relatively large differences in lichen species composition between two sampling periods, but this shift took place along both first DCA axes. Despite overall large compositional dissimilarities at the level of tree phorophytes over time, shifts in lichen species composition were rather phorophyte-specific. We found that in the resampled dataset, the composition of lichens growing mainly on oak and ash was characterized by higher values of vectors representing FRic, species richness, Shannon index, and the CWM of ascomata area. We also found that resampled lichen species associations growing on maple, hornbeam, hazel, and lime were characterized by higher values of FDiv, FDis, and CWMs of WEIV.T and ascospore volume (Figure 4; Table S6). Regarding pine and spruce, the species composition of lichen biota was very similar in both study periods.
Considering the results of the detrended correspondence analysis performed for the substrate level, we found no directional changes in lichen species composition over time, indicating their high substrate-specificity. The only substrates for which we demonstrated relatively clear shifts in lichen species composition were 'branches of living trees' and 'bark of tree trunks or branches of fallen trees.' We found that the resampled lichen species composition of these substrates was characterized by high values of vectors representing the CWMs of light and nitrogen WEIVs ( Figure 5; Table S6). Forests 2021, 12, x FOR PEER REVIEW 12 of 22 Figure 3. Comparison of the changes in parameters of functional diversity (functional richness, evenness, dispersion, and divergence) and CWMs (community-weighted means) of lichen functional traits (ascomata area, ascospore shape, and ascospore volume), CWMs of Wirth's ecological indicator values (light, temperature, continentality, moisture, and nitrogen) and parameters of taxonomical diversity (species richness and Shannon index) between the two periods of lichen sampling (1987-1989 and 2014-2015) were obtained for coniferous, mixed deciduous, and wet forests using linear mixed effect regression. The box covers the 95% confidence range. The thick horizontal line is the median. Black points are outliers. Significance of results: *** p < 0.001, ** p < 0.01, * p < 0.05, n.s.-not significant. Considering the results of the detrended correspondence analysis performed for the substrate level, we found no directional changes in lichen species composition over time, indicating their high substrate-specificity. The only substrates for which we demonstrated relatively clear shifts in lichen species composition were 'branches of living trees' and 'bark of tree trunks or branches of fallen trees.' We found that the resampled lichen species composition of these substrates was characterized by high values of vectors representing the CWMs of light and nitrogen WEIVs ( Figure 5; Table S6).  Table  S1). Arrows represent lichen biota composition characteristics significantly (p < 0.05) correlated with the ordinations results: FRic-functional richness; FDiv-functional divergence; FDis-functional dispersion; AscArea-CWM (community-weighted mean) of ascomata area; SporVol-CWM of ascospore volume; WEIV.L-CWM of light WEIV (Wirth ecological indicator value); WEIV.T-CWM of temperature WEIV; WEIV.C-CWM of continentality WEIV; Rich-species richness; Shan-Shannon index.

Discussion
Lichen species' optimum values changed in Białowieża Forest over a 30-year period for different environmental gradients. Of the lichens for which species' optimum values shifted, the most were those that recently co-occurred with species of higher indicator values for nitrogen. This group included lichens typical of all investigated forest communities, e.g., Zwackia viridis (Ach.) Poetsch & Schied. from mixed deciduous forests, Lobaria pulmonaria (L.) Hoffm. from wet forests, and Calicium salicinum Pers. from coniferous forests. Increased nitrogen deposition in Białowieża Forest has probably contributed to the eutrophication of habitats [26,27], leading to the biotic homogenization of forest communities [41,66]. It may also have caused the transformation of lichen biota towards the increasing proportion of species characterized by higher nutrient requirements. The larger proportion of nitrogen-demanding species was confirmed in all studied forest types, with 43%, 26%, and 23% of variance explained over time in wet, mixed deciduous, and coniferous forests, respectively. However, 8%, 26%, and 27% of variance, in wet, mixed deciduous, and coniferous forests, respectively, were explained by random effects, suggesting the spread of some rare species, recorded in the past only on a few sites, benefiting from potential habitat changes in the investigated area. Such species occurred in adjacent forest types at lower frequencies, e.g., Coenogonium pineti (Schrad. ex

Discussion
Lichen species' optimum values changed in Białowieża Forest over a 30-year period for different environmental gradients. Of the lichens for which species' optimum values shifted, the most were those that recently co-occurred with species of higher indicator values for nitrogen. This group included lichens typical of all investigated forest communities, e.g., Zwackia viridis (Ach.) Poetsch & Schied. from mixed deciduous forests, Lobaria pulmonaria (L.) Hoffm. from wet forests, and Calicium salicinum Pers. from coniferous forests. Increased nitrogen deposition in Białowieża Forest has probably contributed to the eutrophication of habitats [26,27], leading to the biotic homogenization of forest communities [41,66]. It may also have caused the transformation of lichen biota towards the increasing proportion of species characterized by higher nutrient requirements. The larger proportion of nitrogendemanding species was confirmed in all studied forest types, with 43%, 26%, and 23% of variance explained over time in wet, mixed deciduous, and coniferous forests, respectively. However, 8%, 26%, and 27% of variance, in wet, mixed deciduous, and coniferous forests, respectively, were explained by random effects, suggesting the spread of some rare species, recorded in the past only on a few sites, benefiting from potential habitat changes in the investigated area. Such species occurred in adjacent forest types at lower frequencies, e.g., Coenogonium pineti (Schrad. ex Ach.) Lücking & Lumbsch, Melanohalea exasperatula (Nyl.) O. Blanco, A. Crespo, Divakar, Essl., D. Hawksw. & Lumbsch, and Physcia adscendens (Fr.) H. Olivier, but they could spread onto new habitats in subsequent years. These results could be supported by the currently higher functional richness in wet forests, mixed deciduous forests, and coniferous forests, with 74%, 53%, and 25% of variance explained by time and 7%, 2%, and 50% of variance explained by random effects, respectively. Thus, the relatively higher importance of random effects in coniferous forests linked, for instance, with different forest vegetation types in the surroundings, may be evidence that species of higher nutrient demands could shift from neighboring plots and successively spread. Therefore, the current higher functional richness may be positively correlated with shifts in biota towards a larger proportion of more nitrogen-demanding species. Our findings concerning the impact of eutrophication on lichen biota were consistent with results obtained for epiphytic lichen biota [14] and herb layer vegetation [29,41]. In forested mountains, an increased N deposition promotes eutrophic lichens of higher nitrogen preferences, which become dominant over oligotrophic lichens connected with nutrientpoor habitats [20]. The results of the same study also indicated that lichen communities occurring in the driest habitats are more nitrogen-demanding than those occurring in wetter climates at the same N air pollution levels. The shift in the lichen communities on trees twigs in different forest types towards higher contributions of nitrogen-demanding species was reported by Wolseley et al. [67] and Vilsholm et al. [68]. Increasing nitrogen pollution drives gain in the proportion of lichens with specific functional traits, e.g., lichens with chlorococcoid photobiont, foliose species with narrow lobes, and those producing soredia as the main reproduction diaspores [8]. Some authors have employed functional traits of lichens as indicators for assessing environmental changes. Among such traits, crustose, fruticose, and foliose thalli have been mentioned, as has preference to substrate fertility, which can be used to assess transformations in an urban environment [69]. In the case of the examination of conditions in the forest environment, types of growth forms (crustose, narrowly lobed and broad-lobed foliose) and, to lesser extent, reproductive strategy (sorediate, isidiate, or sexual reproduction) functional traits have been mentioned to be the best indicators [19].
Some species co-occurred with species of higher indicator values for temperature. Most of them were lichens typical of mixed deciduous forests growing on hornbeam, ash, oak, and lime, e.g., Ramalina pollinaria (Westr.) Ach., Pseudoschismatomma rufescens (Pers.) Ertz & Tehler, and Lecanora carpinea (L.) Vain., as well as in lesser numbers species typical of coniferous forests, e.g., Chaenotheca chrysocephala (Ach.) Th. Fr. and Micarea melaena (Nyl.) Hedl. However, when considering changes of the increasing proportion of species with higher temperature requirements at the community level, these trends were weak and not significant. In a study on climate warming influence on epiphytic lichen biota, these changes were also not captured [14]; however, only epiphytes growing on trees were analyzed, and epixylic and epigeic lichens were excluded. This could have caused the insignificance of the observed changes due to the limited amount of analyzed data. The analyses performed in the study by Łubek et al. [14] were also of different character: they were based mainly on lichen species compositions and their frequency, contrary to the present analyses based on lichens functional traits and species co-occurrences shifts alongside abiotic (expressed by WEIVs) and biotic (expressed by functional traits) environmental gradients. The present results obtained at the species level indicated that some lichens of higher temperature requirements probably respond to rising temperature values by increasing their frequency. This suggests that while investigating changes in lichen species composition under anthropogenic factors, such as climate warming, the scale of individual species seems to be more appropriate to detect potential effects over time. Shifts in lichen biota considered at larger spatial scales, e.g., at the forest community or tree phorophyte scales, can be overwhelmed by numerous random factors like high site-specific variability or high species turn-over rates [70,71].
We identified changes in the realized optimum occurrence of species co-existing with species of lower indicator values for the continentality environmental factor, which corresponded to the temporal shifts revealed at the species level for the temperature gradient. Among these lichens, most were species characteristic for coniferous forests, e.g., Chaenotheca ferruginea (Turner) Mig., Micarea melaena, and Hypocenomyce scalaris (Ach.) M. Choisy. However, an increase in the proportion of species of low continentality demands was not confirmed at the forest community level over the sampling period-it was instead explained by random effects. This was probably the result of the expansion of hornbeam and lime in Białowieża Forest caused by changes in the environment, including climate change, increasing N deposition, ash dieback, spruce decline, and ungulate herbivory activity, overall resulting in the homogenization of forest communities [41,66]. The colonization of coniferous forests by these broadleaved species, mainly in places adjacent to mixed deciduous forests, may cause shifts in the lichen biota of the coniferous forests towards a greater proportion of species with lower requirements to the continentality environmental factor, which is more typical for deciduous forests.
For several lichens, we recorded changes in the realized optimum for co-occurring with species of lower moisture requirements. Despite the fact that only a few species demonstrated significant shifts in their co-occurrence over time, we also confirmed similar patterns at the community level for all forest types, though mainly for the coniferous (dry) and wet (hydrogenic) forests. This transition in biota was probably caused by the current lower precipitation combined with higher temperatures [26,27], shorter periods of snow cover [28], decreasing groundwater levels [72], and the resultant habitat drying. Our previous paper [14], exclusively concentrating on epiphytic lichens biota during the same time interval, did not detect any significant shift in this respect. However, the present analyses of the entire biota of epiphytic, epixylic, and epigeic lichens showed pronounced, but still relatively weak, rates of changes related to the decreasing humidity of habitats in coniferous and wet forests, as expressed by a low amount of total variance explained by time and random effects. This may suggest that shifts in the moisture optimum for some lichens towards lower values may express initial stages of habitat drying, corresponding to an overall precipitation decrease accompanied by increases in temperatures in Białowieża Forest [26,27,40] and a decrease in snow cover duration [28]. This underlines the importance of studying the full diversity of selected groups of organisms, as they provide more reliable results. Smaller changes towards a greater proportion of lichens with lower moisture requirements were recorded in the mixed deciduous forests. The specific microhabitat with certain humidity and temperature conditions occurring in mixed deciduous forests in Białowieża Forest maintains a relatively stable microclimate, which is characteristic of primeval forests with a long ecological continuity [73][74][75]. This microclimate, to some extent, buffers the effects of the global warming [75], which was confirmed in previous analyses on epiphytic lichens [14], but it does not stop the modifications caused by increased nitrogen pollution input and the resulting progressive eutrophication.
Some species have changed their optimum towards co-occurrence with species of higher values of light requirements, although this trend was only weakly explained by time and random effect in wet forests. The shift in the lichen biota of wet forests is probably associated with the dying of their dominant tree species, ash [76] and alder [77], which transforms forest structure and contributes to greater site overexposure [78], as confirmed by the higher proportion of light-demanding species of epiphytic lichen biota in wet forests [7]. This may also be explained by the close proximity of more open coniferous communities, due to spruce dieback [79], in which light-demanding species of wide ecological scales develop and spread to wet forests due to change in habitat conditions.
We found that some species shifted their optimum in co-occurrence with species with a larger ascospore volume as a functional trait. However, this result was weakly confirmed at the community level only in coniferous and wet forests. Among lichens changing their realized optimum, most species were characteristic for coniferous forests, e.g., Micarea melaena and Lecanora conizaeoides Nyl ex. Cromb. It is likely that species with a larger ascospore volume, which now occur in coniferous forests, have spread from mixed deciduous forests (compare [17]). This can be explained by the appearance of hornbeam tress in some places in coniferous communities bordering with mixed deciduous forests. On these trees, lichen species with a larger ascospore volume typical for hornbeam appeared, e.g., Graphis scripta (L.) Ach., Lecanora pulicaris (Pers.) Ach., and Diarthonis spadicea (Leight.) Frich, Ertz, Coppins & P. F. Cannon.
We found a contrasting pattern for some species characteristic, mainly for wet forests, e.g., Menegazzia terebrata (Hoffm.) A. Massal. and Felipes leucopellaeus (Ach.) Frisch & G. Thor, and coniferous forests, e.g., Placynthiella icmalea (Ach.) Coppins & P. James and Chaenotheca furfuracea (L.) Tibell, which now co-occur with lichens of a lower ascospore shape index. These changes were weakly confirmed at the community level over time. Mass ash dieback [7,76], as well as alder [77] and spruce decline [79,80], in different forest communities created new microhabitats in highly insolated gaps, e.g., lying dead trees, wood of logs, and stumps, which could be rapidly colonized by species characterized by a low ascospore shape index and adapted to a faster spread [15,17,81]. In Białowieża Forest, lichens with a lower ascospore shape index are mainly associated with coniferous communities characterized by harsher microclimatic conditions, such as large temperature and humidity fluctuations, in which habitat filtering is the main ecological mechanism that shapes lichen communities [17]. However, when analyzing the proportion of lichens of the ascospore shape functional trait in the past and at present, we noted a clear change and a greater proportion of species with smaller ascospores. This could be associated with the increasing disturbances in forest stands and resulting decreasing humidity. It is possible that species with a lower ascospore shape index benefit from changes in the environment, are more competitive, and achieve easier reproductive success [81].
It is more difficult to explain shifts in lichens currently showing co-occurrence with species characterized by a larger ascomata area functional trait, as it was not confirmed at the forest community level. These lichens were characteristic for wet forests or mixed deciduous forests, e.g., Pertusaria flavida (DC.) J. R. Laundon, Reichlingia leoppoldi Diederich & Scheid., and Peltigera praetextata (Sommerf.) Zopf. Perhaps this is related to forest dynamic processes, like long tree natural rotation and the close placement of newly regenerated trees forming greater microhabitat diversity, which affected the spreading success of some lichens with a larger ascomata area [71]. This might also be related to the achievement of greater internal microclimatic stability in terms of temperature and humidity, which was found for mixed deciduous forests in the study area [14,17]. Further research is needed to explain the obtained results.
In analyzing changes in lichen biota at the phorophyte level, within different parameters of functional diversity, we observed only significant increase of functional richness. This took place on different phorophytes, especially on ash and oak, which are known in the study area to have the highest lichen species diversity both in the past and at present [7,45]. The increase in the functional richness on these trees was only correlated with the increased proportion of lichens with a larger ascomata area. On the phorophyte level, opposite to the community level, we did not observe changes in the proportion of nitrogen-demanding species or in the contribution of less moisture-demanding species.
We did not record any clear changes in the functional richness of lichen biota at the substrate level. This confirms the significant specificity of lichens to particular substrate types [82][83][84]. It also may indicate relatively stable conditions in microhabitats forming within certain substrates, e.g., the foot of trees, bark of trees, and root system of fallen spruces. Small changes toward the higher proportion of nitrogen-demanding species were only recorded for the branches of living trees, the bark of tree trunks, and the branches of fallen trees, which was in line with the findings of our previous study [14]. This could have been a result of greater accessibility of nitrogen compounds in the higher parts of the tree crown due to bird droppings (e.g., [69,85]). However, certain works have pointed out that epiphytes growing on trees can be used as an early warning system to detect nitrogen deposition [67,68]. Our results may also have been biased by potentially higher sampling efforts in cases of lichens growing on these types of substrates during our investigations.

Conclusions
Changes in the lichen biota in Białowieża Forest were evident at the level of individual species, such as shift of species' realized optimum. These changes were not always detectable in the forest community, phorophyte, and substrate levels, though they were clearly indicated in biota's shifts in the case of a higher proportion of nitrogen-demanding species, and they were relatively weak in the lower proportion of moisture-demanding lichens. Though we did not have direct climatic data over the time period, our results suggest that the modification of lichen biota took place due to changes in the environment under the impact of climate warming and eutrophication. These transformations are probably in the initial stage, and that is why they are not clearly detectable at the level of forest communities. These transformations may also be related to species turn-over and the site-specificity of microclimates, and they thus may demonstrate the high resistance of the primeval forest ecosystem to indirect human impact. We also suppose that other overlapping factors may have contributed to the observed shift in biota composition, e.g., a decrease in the abundance of ash and spruce and an expansion of lime and hornbeam. Nevertheless, in Białowieża Forest, the influence of climate change on lichens can be tracked down.
Analyses of the shifts in the species composition of lichen biota provide generalized information on the rates of these changes and the importance of anthropogenic factors. Studying the time dynamics of the lichen species' realized optima can be used for the indirect identification of shifts in lichen biota as a response to changes in the forest environment. Concurrent analyses of the changes of species optimum and functional diversity, contrary to simple analyses of shifts in species composition alone, can better reveal the directions of changes in lichen biota over time.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12060686/s1. Table S1: List of species recorded in the study area.; Table S2: Parameters of variables passively fitted to results of the PCA ordinations showing temporal shifts in composition of lichen species at forest community level. Determination coefficients R 2 and p-values were estimated using permutation tests with 999 iterations. Abbreviations: CWM-community weighted mean; WEIV-Wirth ecological indicator value. Significant results are in bold; Table S3: Parameters of mixed effect models showing temporal shifts in lichen species' functional diversity components, community weighted means of functional traits, Shannon index and species richness regarding coniferous forest communities. For testing shifts in species richness over time, the generalized mixed effects model with Poisson distribution was employed, while regarding remaining parameters of lichen species composition the linear mixed models were used. Abbreviations: CWM-community weighted mean; WEIV-Wirth ecological indicator value. Significant results are in bold; Table S4: Parameters of mixed effect models showing temporal shifts in lichen species' functional diversity components, community weighted means of functional traits, Shannon index and species richness regarding mixed deciduous forest communities. For testing shifts in species richness over time, the generalized mixed effects model with Poisson distribution was employed, while regarding remaining parameters of lichen species composition the linear mixed models were used. Note that regarding functional dispersion, the model assumptions were not met due to the very small variance of random effects. Abbreviations: CWM-community weighted mean; WEIV-Wirth ecological indicator value. Significant results are in bold; Table S5: Parameters of mixed effect models showing temporal shifts in lichen species' functional diversity components, community weighted means of functional traits, Shannon index and species richness regarding wet forest communities. For testing shifts in species richness over time, the generalized mixed effects model with Poisson distribution was employed, while regarding remaining parameters of lichen species composition the linear mixed models were used. Note that regarding functional evenness, functional divergence and spore volume, the model assumptions were not met due to the very small variance of random effects. Abbreviations: CWM-community weighted mean; WEIV-Wirth ecological indicator value. Significant results are in bold; Table S6: Parameters of variables passively fitted to results of the DCA ordinations showing temporal shifts in composition of lichen species at tree phorophyte and substrate levels. Determination coefficients R 2 and p-values were estimated using permutation tests with 999 iterations. Note the presence-absence data transformation prior to DCAs. Abbreviations: CWM-community weighted mean; WEIV-Wirth ecological indicator value. Significant results are in bold.
Author Contributions: Conceptualization and data collection, A.Ł. and M.K.; software, data curation, and analysis P.C. and A.Ł.; writing-original draft preparation, A.Ł., M.K., B.J., and P.C.; supervision and project administration, B.J. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare that they have no conflict of interests.