A biogeographical study of red listed lichen species at temporal and spatial scales within protected and non-protected areas

The present study is focused on the temporal and spatial distribution of red listed lichen (RLL) species identified in both non-protected areas (NPAs) and protected areas (PAs) in Romania. This study revealed different scenarios of RLL based on two major patterns: (1) the fate of RLL before and after their designation in the red list in Romania and (2) the fate of RLL before and after the designation of PAs in Romania. Generally, the occurrences of RLL were well represented in time and space in both NPAs and PAs through geomorphological and biogeographical units. In particular, the occurrences of RLL were well represented, especially in hilly areas within PAs before and after their official designation, and this was an important aspect over a long period of time. Although NPAs were not less important regarding the occurrences of RLL species, it was observed that they followed the same pattern as PAs, namely, they were well represented in hilly areas and over a long period of time. The bioregions were significant for RLL species in both NPAs (continental and stepic) and PAs (continental, panonic, and stepic). As a significant finding of this study, NPAs and PAs should be subjected to an adequate conservation regime due to their biotical traditional heritage.


Results
The RLL species were identified in 92 sites within PAs and in 99 sites within NPAs. Among a total of 18 RLL species, all were found in the PAs, while in the NPAs, 15 RLL species were identified (supplementary file: Table S1 and Table S2). The distribution and number of RLL species within the geomorphological units and biogeographical regions in both NPAs and PAs are given in supplementary file: Table S3. The occurrence of each RLL species in the geomorphological unit and biogeographical region in the four periods of time are presented within the NPAs (supplementary file: Table S4) and PAs (supplementary file: Table S5).
Regarding the statistical analysis, it was important to control confounders based on one-way NPMANOVA for RLL species data within NPAs and PAs. After removal, the data belonging to confounders such as Bryoria lanestris (Ach.) Brodo & D. Hawksw. and Nephromopsis chlorophylla (Willd.) Divakar, A. Crespo & Lumbsch, the results of this test based on the adjustment of confounders did not indicate significant associations between datasets represented by the other RLL species within the NPAs and PAs (F = 3.76; p = 0.06). This indicates that the results from the following statistical analyses were not biased by associations between the datasets.
In the present study, the RLL species number did not significantly differ among the NPAs and PAs (chi 2 = 2.75; p = 0.99).
Within the NPAs, the time periods 1850-1900 and 2001-2020 were excluded from statistical analyses due to spatial autocorrelation results, because of these time periods and their associated geomorphological and biogeographical units were biased by associations between the datasets. Also, within the PAs, the time periods 1850-1900 and 1901-1950 were excluded from statistical analyses based on the same statement as above.
In the seriation analysis, the RLL species were ordered across geomorphological and biogeographical units throughout the time periods. Thus, throughout the time periods, the ordered pattern of RLL species across the geomorphological and biogeographical units was significant within both the NPAs and PAs (Table 1). Within the NPAs, the RLL species were ordered in the hilly areas within continental (CON) during the second period of time and within stepic (STE) and CON bioregions during the third period of time (Table 1). Within the PAs, the RLL species were ordered in the third and fourth periods of time in the hilly areas situated in the CON and panonic (PAN) bioregions (in the third periods of time) and in the STE and CON bioregions in the fourth period of time (Table 1).
The occurrences of RLL species changed along spatial and time attributes. Thus, in the NPAs, major changes over time periods (1901-1950 and 1951-2000) in the occurrences of RLL species were found in geomorphological units such as hills and their corresponding CON and STE biogeographical regions (Table 2). In the NPAs, the  Table 2. The ordering of RLL species in the geomorphological units and biogeographical regions over time within NPAs (according to the seriation analysis performed in PAST, based on the unconstrained algorithm). 2-1901-1950; 3-1951-2000; asterisks indicate the occurrence of RLL species in the geomorphological units and biogeographical regions over time; empty cells indicate no available data.

STE 3 Hill CON 2
Ramalina obtusata (Arnold) Bitter * * * * RLL species were concentrated mostly in the CON bioregion during the second half of the twentieth century ( Table 2). In the PAs, main changes over time periods (1951-2000 and 2001-2020) in the occurrences of RLL species were found in geomorphological units, such as hills and their corresponding CON, PAN and STE biogeographical regions ( Table 3).
The Kruskal-Wallis test revealed significant differences between occurrences of RLL species for NPAs within geomorphological units and biogeographical regions over time periods (chi 2 = 20.64; p = 0.0001) and for geomorphological units and biogeographical regions within PAs (chi 2 = 15.55; p = 0.001).
Post hoc comparison tests pointed out significant differences in RLL species occurrences across biogeographical regions during the second and the third periods of time (STE versus CON, p = 0.0008) for NPAs. Within hilly areas from NPAs, significant differences in RLL species identified in the STE bioregion were recorded during the second and the third periods of time (hill versus STE, p = 0.0007).
Post hoc comparison tests performed for PAs indicated significant differences between pairwise data represented by RLL species across both geomorphological units and biogeographical regions during the third and the fourth periods of time. Furthermore, in the PAs, significant differences were obtained based on post hoc Mann-Whitney pairwise comparisons of RLL species occurrences within biogeographical regions during the third and the fourth periods of time (STE versus CON p = 0.003; CON versus PAN p = 0.0008). During the third and the fourth periods of time within hilly areas from PAs were recorded significant differences in RLL species identified in the STE and PAN bioregions (hill versus STE, p = 0.003; hill versus PAN, p = 0.0007).
In the NPAs, the occurrences of the RLL species within the geomorphological units and biogeographical regions were significantly different over the time periods but not in the periods of time when were designated in the Red List of lichens in Romania (supplementary file: Table S6). Otherwise, regarding the impact of the designation of PAs on RLL species, it was demonstrated that in PAs, the RLL species showed significant differences during the third  and fourth (2001-2020) periods of time within the geomorphological units and the biogeographical regions. Consequently, the RLL species showed significant differences in their occurrences before and after legal designation of the PAs (between 2000 and 2010) in Romania. In the time period 2007-2013, the RLL species for Romania were published and designated (supplementary file: Table S6). Based on the information presented in supplementary file: Table S6, 18 RLL species were designated in Romania 21,22 . In the PAs, the RLL species were significantly differentiated within geographical attributes before and after their designation in the two Red Lists (Tables 1, 3 One-way ANOSIM indicated significant differences in the RLL species between the NPAs and PAs (R = 0.41; p = 0.006) based on their occurrences in the geomorphological units and biogeographical regions over the various time periods (Fig. 1).
The post hoc pairwise ANOSIM test indicated significant differences (p = 0.006) in the occurrences of RLL species between all pairs of the two studied groups defined by the NPAs and PAs in their geomorphological units and biogeographical regions over the time periods.
The statistical analysis based on the Kolmogorov-Smirnov test indicated nonsignificant differential collection rates during all periods of time and their associated biogeographical and geomorphological variables within the NPAs and PAs. The one-way NPMANOVA did not indicate significant differential collection rates between the NPAs and PAs during all periods of time and their associated biogeographical and geomorphological variables.

Discussion
In the present study the species number among NPAs and PAs was not significantly different. However, other studies showed that the species number was better represented in PAs than in NPAs, due to effective control of human activities 5 . Nevertheless, NPAs have particular importance for biodiversity conservation and support measures for sustainable management and limitations of anthropogenic pressure should be implemented 2,4,14 .
Usually, red lists include threatened and extinct species with great importance for environmental resource conservation 4,26,27 . Red-listed lichen species are related to ancient habitats that contain specific microhabitats 28,29 and therefore it is necessary to secure RLL survival by adequately protecting these habitats 10 . The occurrence of Table 3. The ordering of RLL species in the geomorphological units and biogeographical regions over time within the PAs (according to the seriation analysis performed with the PAST software, based on unconstrained algorithm). 3-1951-2000; 4-2001-2020; asterisks indicate the occurrence of RLL species in the geomorphological units and biogeographical regions over time; empty cells indicate no available data.  34 . At small spatial scales, geomorphological units such as hilly areas significantly affect the variation in the occurrences of RLL species; therefore, geomorphology associated with different climatic conditions (mountain habitats with low temperature and high amounts of precipitation are different from hilly and plain habitats with high temperatures and low amounts of precipitation) is an important driver of various tree species 30 that harbour RLL species due to tree diversity and microhabitat quality 3,10 . At a large spatial scale, geomorphology and biogeographical regions and their associated abiotic parameters play an important role in determining the occurrence of lichen species across various habitats 13 . It is worth considering that even at higher altitudes in mountain areas, the anthropogenic pressure related to a weak efficacy of legislative systems has a negative impact on natural habitats 17 , which leads to many difficulties in biodiversity conservation 31 . PAs should be represented by a wide range of geomorphological features and the sustainable connectivity needed for the long-term conservation of biodiversity 32 .
The observed distribution of RLL species occurrences across landscape features, such as the structural heterogeneity of traditional habitats in NPAs and PAs, is attributed to the lack of anthropogenic impact 33 . Within these remote areas, the main management land use type is based on traditional practices that do not represent a major negative impact on the environment 33 . In Romania, many NPAs are situated in remote sites with long-term continuity and connectivity of different types of habitats, such as old-growth forests, meadows, rocks, wetlands, etc. 34 . The structural complexity of the traditional landscape represents suitable habitats for lichen species and therefore should be a focus for the establishment of PAs in Romania 35 . Habitat connectivity among protected landscapes and their heterogeneity are important criteria, especially for threatened species and their habitats 36,37 .
Within this study, the RLL species were significantly represented during various periods of time across NPAs and PAs. As in this study, both NPAs and PAs are important for RLL species 6 , which emphasizes their importance in biodiversity conservation 38 . The different geographical attributes (geomorphology and biogeographical regions) that characterize PAs support a high richness of lichen species 8 . In this regard, the territory of Romania is well represented by geomorphological units 30 and biogeographical regions 34 , and RLL species have been identified in various vegetation communities 22,[39][40][41][42][43] .
In many regions of the world, ancient habitats are restricted and isolated areas 8,44 surrounded by a matrix of strongly disturbed anthropogenic areas 27,38 , which significantly affects lichen communities 43,45 . The connectivity of habitats is an important attribute for biodiversity conservation within PAs 44 during dynamic landscape stages and it impacts species assemblages and the quality of their habitats 46 .
The RLL species represent a natural heritage because some of them, such as Lobaria pulmonaria (L.) Hoffm., are widely recognized as valuable indicators of ancient forests threatened by industrial activities and the management of forested areas 47 . Furthermore, L. pulmonaria is an RLL species in many countries in Europe 47,48 , and it is listed as endangered in Europe 49 .
The support of biodiversity conservation is highlighted by a complex frame of concepts and strategies based on historical arguments with their origins in scientific and political resolutions 50 . It is important to grasp that the designation of PAs does not represent a universal remedy for biodiversity conservation; consequently, biodiversity conservation within PAs depends on sufficient ecological resources secured by low human pressure 27 . In European countries, the conservation of biodiversity in the current context could not completely eliminate www.nature.com/scientificreports/ anthropogenic impacts, but it could be reduced so that the species and habitats are not affected or lost 51 . The common efforts of politicians and local communities to advance strategies for long-term environmental sustainability are key to the success of nature conservation efforts 50,51 . Conclusions. Non-protected areas are important for RLL species; therefore, these areas should be subjected to adequate management for the long-term conservation of habitats and their associated RLL species. Additionally, PAs and their geographical attributes have an important role for RLL species. Through time, the occurrence of RLL species has been significantly changed, and particular geographical patterns have emerged that support the habitat of these important species in NPAs and PAs. The habitats represented by NPAs and PAs are both equally important for conservation of the RLL species in Romania. The time periods clearly suggest that it does not matter when the PAs were designated but how well their traditional natural heritage was conserved. This statement is supported by the fact that the RLL species survived over time, and more importantly, these species are still found within NPAs and PAs in Romania.

Methods
Legislative framework. The RLL species from Romania are not subjected to any legislative framework.
The designation of PAs is mainly enforced by Law no. 5/2000 23 as a tool of national legislative support for habitats and species conservation. Additionally, the other PAs were designed by the Natura 2000 Network 24,25 .
Data sources. Data on the time and spatial occurrences of RLL species were obtained from the Mycological Herbarium of the Institute of Biology from Bucharest (BUCM), which belongs to the Romanian Academy (supplementary file: Table S1 and Table S2), and the author's unpublished and published data (supplementary file: Table S1 and Table S2). For this study, all RLL species from Romania were considered, for a total of 18 species (supplementary file: Table S6). The nomenclature of the lichen species is according to www. index fungo rum. org. Geographical distribution of RLL species in Romania. The RLL species were georeferenced based on a geographical information system represented by geographical coordinates measured in decimal degrees within the studied NPAs and PAs, which were widely spread across a range of different geomorfological units, such as plains, hills, and mountains. The old chronological data included in supplementary file: Table S1 and Table S2 of RLL species were represented by toponyms and localities (towns and villages). The geographical distribution of RLL species used in this study in the case of old chronological data, which includes toponyms and localities, was based on the geographical information system represented by the Universal Transverse Mercator (UTM) grid system of 100 × 100 km for Romania 52,53 . A universal transverse mercator grid was used to identify toponyms and localities in Romania in the case of old chronological data; most toponyms and localities were published by 52,53 based on a UTM grid. The recent chorological data included in supplementary file: S1 and Table S2 taken into account in this study were also represented by toponyms and localities and were geographically referenced using GPS coordinates based on degrees, minutes, and seconds. An exception is the newer records that were determined using a GPS set to geographical coordinates measured in decimal degree to give a more precise location for plotting. The spatial resolution of some older data represented by toponyms is 100 km square whilst the spatial resolution of newer data is represented by a few meters due to the accuracy of handheld GPS. In the database of this study, the entire dataset represented by old and recent chronological data was converted from UTM geographical coordinates and geographical coordinates given as degrees, minutes, and seconds in the decimal degree World Geodetic System of 1984 (WGS84) using the option Military Grid Reference System (MGRS)/UTM reference (UTMREF) based on WGS84 available at the following site https:// coord inates-conve rter. com.
A map that shows the geomorphological units and biogeographical regions from Romania (Fig. 2) was generated based on the decimal degree (WGS84) geographical coordinates in ArcGIS 10.4 54 using the ETRS89 LAEA reference system. Elevation is represented based on the 30 m digital elevation model over Europe freely available from the Copernicus Land Monitoring Service website (https:// land. coper nicus. eu). Biogeographical regions were drawn using the shapefile available on the website of the Romanian Ministry of Environment, Water and Forestry (http:// www. mmediu. ro) using the colour codes specified by the European Environmental Agency on their website (https:// www. eea. europa. eu).
Within this work, in the case of old records (published data, especially in the last two centuries), it is not known how the data were collected by authors to avoid spatial biases such as data collected near roads due to www.nature.com/scientificreports/ easy accessibility. Based on this idea, in the author's original database, most data on the geographical distribution of RLL species in Romania were collected in remote areas (far from roads and localities) to avoid spatial accessibility biases. There were some exceptions in the case of one specimen of Hypotrachyna sinuosa (Sm.) Hale, which was collected near a national road in the southern part of Romania (all of the other specimens of H. sinuosa considered in this study were collected in remote areas). Additionally, another exception was the case of three specimens of Cetraria islandica subsp. islandica (L.) Ach. collected near a national road in the central part of Romania. In regards to the old records, because it is not known how the data were collected, to avoid spatial biases (near roads and certain localities) a statistical analysis in this sense was impossible.
Statistical analysis. Lichen occurrence assessment was established as the presence/absence of a species recorded within geomorphological units and biogeographical regions through time periods within PAs and NPAs.
In the statistical analysis, an adjustment was performed to control all of the confounders so that the results of the statistical analyses reflected true statements 55 . In this regard, a one-way NPMANOVA test was performed to control for the confounders 55 . In the one-way NPMANOVA test, the distance measure was based on the Kulczynski index appropriate for binary data. The test was computed by 9999 permutations, and finally, a post hoc test was provided between the pairwise of two groups represented by PAs and NPAs 56 . The one-way NPMANOVA test identified two confounders based on the data of two RLL species, Bryoria lanestris (Ach.) Brodo & D. Hawksw. and Nephromopsis chlorophylla (Willd.) Divakar, A. Crespo & Lumbsch, and as a consequence, the data for these species were eliminated from all statistical analyses.
Spatial autocorrelation was used to identify the dependence among the spatial data (geographical coordinates of the RLL species) and the biogeographical and geomorphological data of the RLL species within both NPAs and PAs 57 . Because the dataset is represented by binary data (presence/absence of RLL species), it is recommended to analyse spatial autocorrelation by using generalized estimating equations 58 . The generalized estimating equations were performed with R software 59 by using the geepack package [60][61][62] . The generalized estimating equations indicated a significant correlation within the spatial data represented by geographical coordinates and data recorded for plain, mountain, and alpine variables within NPAs and PAs; therefore, all of these spatial data and their associated RLL species were removed from Seriation, Kruskal-Wallis, one-way ANOSIM, and one-way NPMANOVA analyses.
Comparisons of the total species number between the NPAs and PAs was performed using the chi-squared test 56 . www.nature.com/scientificreports/ Seriation analysis was used based on a binary data matrix that represents the presence/absence of RLL species within geomorphological and biogeographical units in Romania for both NPAs and PAs. Seriation analysis was applied to identify the chronological pattern in RLL species across geomorphological and biogeographical units. Thus, this analysis was used to reveal the (a) ordination and differences of RLL species along spatial (geomorphological and biogeographical units) and time (semicentenary periods of time) gradients, and (b) how the occurrences of RLL species are ordered along these gradients. The time periods are an important informational support over time because represent a time landmark for RLL species occurrences from spatial point of view. Furthermore, seriation analysis was used to highlight the changes in RLL species occurrences along biogeographical and geomorphological gradients throughout the time periods. The seriation analysis was performed using constrained and unconstrained algorithms 56 . Thus, in the constrained algorithm, the RLL species were chronological ordered (Table 1) whilst in the unconstrained algorithm, the RLL species were ordered along the representative period of time and its associated geomorphological and biogeographical units (Tables 2, 3). A Monte Carlo simulation was run using the constrained algorithms 56 .
The Kruskal-Wallis test was used to analyse differences within NPAs and PAs based on all of their attributes, while one-way ANOSIM was used to reveal the differences between PAs and NPAs based on all of their attributes.
The Kruskal-Wallis test was used to identify the differences in the occurrences of RLL species along geomorphological and biogeographical units throughout the time periods within the NPAs and PAs. A post hoc pairwise Mann-Whitney test based on p values was used to compare the occurrences of the RLL species along the spatial and time gradients in the NPAs and PAs 56 .
One-way ANOSIM was used to test the differences between the data represented by NPAs and PAs (studied groups). The significance of the test was computed by permutation of group elements with 9999 replicates. As a distance measure, the Jaccard distance index was chosen since it is appropriate for presence/absence data. A post hoc test was performed to identify differences in all pairs of studied groups represented by NPAs and PAs data 56 .
Different time periods were statistically compared; therefore, collection rates at different time periods were calculated for both NPAs and PAs. The collection rates were calculated as follows: the number of RLL species records in the biogeographical and geomorphological units relative to the unit of time, namely, 50 years 63 . The collection rates for NPAs and PAs were calculated for each period of time taken into account within this study (supplementary file: Tables S7 and Table S8). The statistical analysis of collection rates within the NPAs and PAs was performed based on a distance matrix calculated using the Gower distance index between all pairs of time periods across the biogeographical and geomorphological attributes 56 . The significance of the differential collection rates between all pairs of time periods and their geographical and geomorphological variables was tested using the Kolmogorov-Smirnov test 56 .
One-way NPMANOVA was used to test the differential collection rates between the NPAs and PAs, each of which was represented by a distance matrix that included the periods of time and their biogeographical and geomorphological variables. For each distance matrix, we selected the Gower distance index. Permutation was performed with 9999 replicates 56 .
Statistical analyses such as Seriation, Kruskal-Wallis, one-way ANOSIM, and one-way NPMANOVA were performed with PAST software 56 . Generalized estimating equations were performed with R software 59 .

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.