Testing the Poleotolerance Lichen Response Trait as an Indicator of Anthropic Disturbance in an Urban Environment

Urban environments are densely populated areas buzzing with a wide range of anthropic activities that cause disturbances like air pollution or the heat island effect, threatening both human and environmental health. Mitigating its impacts implies understanding the integrated effects that those disturbances exert on urban environments. Lichen biodiversity is frequently used as an ecological indicator, being able to integrate its effects in a quantifiable way. The poleotolerance response trait classifies lichens according to their tolerance to human disturbance, but it was developed for Italy’s flora and has seldom been applied outside Italy or in urban context studies. The aim of this work was to assess this trait suitability as an indicator of urban anthropic disturbance and test it outside Italy. For that, we sampled lichen diversity in 41 green spaces in Lisbon. Lichens were classified into the respective poleotolerance trait functional groups and their community weighted mean related with three type of environmental variables used as surrogates of urban disturbance. We showed that disturbance-tolerant functional groups could be used as an ecological indicator of the integrated effects of environmental disturbances. Some species were clearly misclassified, so we propose reclassification for those. Natural and seminatural functional groups did not behave as expected. Nevertheless, disturbance-tolerant functional groups have the potential to be used in in other Southern European cities.


Introduction
Anthropic pressures to ecosystems have been rising along with the world's population. According to Roser and Ortiz-Ospina [1], from 1900 to 2015, the world's population increased from 1.5 billion to 6.1 billion. Projections estimate that it can stretch to a maximum of 12 billion by 2050 and 17.3 billion by 2100 [2]. Population growth depletes ecosystems worldwide, as our demands for more resources and services rise. Human activities lead to increased air, soil, and water pollution [3], climate change [4], deforestation [5,6], and others. Combined, these disturbances create a major threat to ecosystems worldwide and all fauna and flora that it comprises.
Coupled with the population growth phenomenon is the increasing amount of people moving from rural to urban areas [7]. According to the UN [8], in 2007, for the first time in the history of humankind, the rural population was surpassed by the urban population. If current trends continue, it is expected that by 2050, 66% of the total world population will be living in urban areas. These mass movements create even more pressure on these already densely populated areas, due to the need for more space for infrastructures, namely housing and industries, public and private transportation, agricultural practices, and others. highly disturbed habitats. This trait was explored in a few other studies in Europe [37,38], one of them conducted in an urban environment [23]. The major potential of this trait is that it does not focus solely on the effects of pollutants in lichens communities, but rather it looks at all anthropic disturbances in an integrated way and also how those disturbances, their type, and intensity have changed over the years. This trait is not intended to distinguish between different human disturbance sources, like other traits (e.g., eutrophication, humidity requirements or bark pH), but as a simple integrating trait that indicates disturbance as a whole. In this sense, it could be a useful quick tool to help to signal overall human disturbance, for instance, signaling areas that deserve attention and where future studies should focus to identify the individual environmental factors driving the impacts. However, it is still unclear if it has the potential to be used in urban environments outside Italy. As Gilbert [39] stated, lichens distribution in urban environments is not exclusively determined by air pollution but by a series of other environmental factors like pH, habitat or climatic conditions.
The main aim of this work was to assess the suitability of the poleotolerance response trait as an indicator of anthropic disturbance in a context of urban environments, particularly outside Italy, for where this classification was designed for. This was tested in the densely populated city of Lisbon, where contrasts between housing areas, green spaces, and industries form a marked gradient of disturbance. This disturbance gradient should be potentially reflected on epiphytic lichen communities in a way that allows us to test this trait efficiency [40,41].

Confirming Species Classification
This study was conducted in 41 green spaces of Lisbon municipality. This municipality is part of the Lisbon Metropolitan Area. This area is inhabited by 2.8 million people [42] and has a population density of 1484/km 2 [43]. The climate is Mediterranean. From 1960 to 2017, Lisbon recorded an annual mean temperature of 17.2 • C and an annual mean precipitation of 779.8 mm [44]. In 2017, the annual mean temperature was of 18.3 • C and the annual mean precipitation was 505.1mm. Approximately 370 thousand cars are estimated to circulate daily in the city [45].
Lisbon green spaces were identified using the European Urban Atlas [46], and more green spaces were added using manual photo interpretation, to ensure green spaces with less than 1 ha were included. The selection of the green spaces to sample followed a stratified sampling design. The stratification was done by location and area. The city was divided into four quadrants to guarantee sampling sites were distributed equally among them. Five classes of size were established (<0.1 ha; 0.1 to 0.5 ha; 0.5 to 10 ha; 10 to 40 ha; >40 ha), and two green spaces were chosen in each, for a balanced distribution. This resulted in 40 sampling sites (4 quadrants × 5 area classes × 2 green spaces). One more site was added to the previous ones. This last site was positioned inside the Parque Florestal de Monsanto (Monsanto Forest Park) due to the fact that this is, by a large extent, the biggest green space in Lisbon (1000 ha). Thus, a total of 41 sites were sampled ( Figure 1).

Lichen Sampling
Macro and micro epiphytic lichens were sampled during March, April, and May of 2015. The four suitable trees closest to the centroid of each sampling site, in each green space, were selected and sampled following the European standard methodology [19]. We opted for using the trees closest to the centroids because they offer the highest potential for disturbance reduction due to the fact that they are the most protected area inside the green space. Although the European methodology [19] recommends using only one species of phorophyte, this was impossible due to the great diversity of phorophytes in Lisbon. We estimate that 60% of phorophytes were Jacaranda mimosifolia D. Don, as this species is one of the most widely spread species in Lisbon green spaces, with the adequate characteristics for this work. Tipuana tipu (Bentham) Kuntze accounted for around 30%, and the rest were from unidentified species, as sampling was conducted during winter, and these remaining phorophytes were all deciduous. To address this limitation, only phorophytes with medium bark roughness were used for sampling.

Lichen Sampling
Macro and micro epiphytic lichens were sampled during March, April, and May of 2015. The four suitable trees closest to the centroid of each sampling site, in each green space, were selected and sampled following the European standard methodology [19]. We opted for using the trees closest to the centroids because they offer the highest potential for disturbance reduction due to the fact that they are the most protected area inside the green space. Although the European methodology [19] recommends using only one species of phorophyte, this was impossible due to the great diversity of phorophytes in Lisbon. We estimate that 60% of phorophytes were Jacaranda mimosifolia D. Don, as this species is one of the most widely spread species in Lisbon green spaces , with the adequate characteristics for this work. Tipuana tipu (Bentham) Kuntze accounted for around 30%, and the rest were from unidentified species, as sampling was conducted during winter, and these remaining phorophytes were all deciduous. To address this limitation, only phorophytes with medium bark roughness were used for sampling.
Sampling was conducted using a grid with 50 cm × 10 cm (divided into five equal squares, 10 cm × 10 cm). The grid was attached at chest height, to the tree trunk facing t he four cardinal points, and the frequency of each present species (maximum of 20 and minimum of 1, per tree) was noted. In addition to bark roughness, phorophyte selection also took into account trunk inclination (no more than ten degrees deviation from vertical), diameter at chest height (between 50 cm and 150 cm), and with no visible trunk damage, knots or other obstacles for epiphytic lichens development. All phorophytes' diameter and coordinates were recorded using a measuring tape and a GPS device, respectively. A magnifying lens was used for in-site species determination. Species that could not be determined in the field were properly collected and stored in envelopes and later identified in a Sampling was conducted using a grid with 50 cm × 10 cm (divided into five equal squares, 10 cm × 10 cm). The grid was attached at chest height, to the tree trunk facing the four cardinal points, and the frequency of each present species (maximum of 20 and minimum of 1, per tree) was noted. In addition to bark roughness, phorophyte selection also took into account trunk inclination (no more than ten degrees deviation from vertical), diameter at chest height (between 50 cm and 150 cm), and with no visible trunk damage, knots or other obstacles for epiphytic lichens development. All phorophytes' diameter and coordinates were recorded using a measuring tape and a GPS device, respectively. A magnifying lens was used for in-site species determination. Species that could not be determined in the field were properly collected and stored in envelopes and later identified in a laboratory resorting to a microscope and different identification keys [47,48]. Species nomenclature follows Nimis [36].
Species were classified according to the maximum value of the poleotolerance response trait available in ITALIC [36] (Table 1). Trait classification and species frequencies were combined to calculate the community-weighted mean (CWM) for each functional group, in each green space. The CWM is a trait-based metric representing a measure of functional structure. As this trait is categorical, the CWM here calculated represents the mean proportion of each functional group in the community, weighted by the abundance of species belonging to those groups [49]. CWM was calculated using the "dbFD" function from the FD package [50] of R CRAN software [51]. Table 1. Description of each of the four functional groups of the poleotolerance response trait. Description follows as in Nimis [36].

Functional Group Description
Natural Species which exclusively occur on old trees in ancient, undisturbed forests. Seminatural Species occurring in natural or seminatural habitats Low disturbance Species occurring in moderately disturbed areas High disturbance Species occurring in heavily disturbed areas

Environmental Variables
Three types of variables were calculated for each green space: (1) The area (m 2 ) of each green space; (2) the normalized difference vegetation index (NDVI), both for the green space centroid (NDVIgs) and for a 100 m buffer around each sampling site (NDVIb100), used as a proxy of vegetation density; and (3) the surrounding land cover area around the green space.
Green spaces areas were determined using ArcGIS Desktop 10.6.1 (ESRI). NDVI values refer to May 2015 (to remove the effects of grasslands, after image correction) and were calculated based on a Landsat 8 image (30 meters resolution). NDVIb100 was calculated so we could take into consideration the possible effects of nearby green spaces. NDVI is used as a proxy of vegetation density and green space area. However, that does not necessarily mean that bigger green spaces have higher NDVI values, as the NDVI values are higher or lower depending on the amount of surface covered by green vegetation and not on the area of the green space (e.g., a bigger green space with trees and bare soil will have lower NDVI values than a smaller green space with trees and grass instead of bear soil). Characterization of the land cover surrounding each sampled green space was also made in ArcGIS Desktop 10.6.1 (ESRI) using the information available in the Urban Atlas. Buffers with seven different radii (50 m; 100 m; 200 m; 500 m; 1000 m; 1500 m; 2000 m) were created around each sampling site. Buffers should be interpreted as drawn imaginary circles, with different radius, created with the sole intention to delimit distinct areas around the sampling sites. These buffers were used to identify all land covers surrounding all the sampling sites, at different distances. A total of 25 land covers were listed in, at least, one of the buffers' distance: Airports; arable land; construction sites; continuous urban fabric (S.L.: >80%); discontinuous dense urban fabric (S.L.: 50-80%); discontinuous medium density urban fabric (S.L.: 30-50%); discontinuous low density urban fabric (S.L.: 10-30%); discontinuous very low density urban fabric (S.L.: <10%); fast transit roads and associated land; forests; green urban areas; herbaceous vegetation associations; industrial, commercial, public, military, and private units; isolated structures; land without current use; mineral extraction and dump sites; open spaces with little or no vegetation; other roads and associated land; pastures; permanent crops; port areas; railways and associated land; sports and leisure facilities; water and wetlands. From the 25 land cover types, five were discarded (airports; construction sites; mineral extraction and dump sites; open spaces with little or no vegetation; and wetlands), as they were not significantly correlated with the CWM nor with the ordination scores, at any buffer distance (Table S1). For the remaining land covers, the best buffer distance was selected. After that, land cover types with the same buffer distance and similar correlation trend (correlation signals) with the first and second axes from the NMS and the CWMs of all three functional groups were grouped. Four groups were created (HDLC; LDLCa; LDLCb; LDLCc).

Data Analysis
Nonmetric multidimensional scaling (NMS) was performed to identify the main drivers of lichens communities' changes. This analysis was run in a relativized matrix of species abundances, to avoid biasing results and impairment of comparisons [16]. Bray-Curtis distance was used, and the best NMS solution was chosen from 50 runs (with real data), each starting randomly (500 iterations per run), and evaluated with a Monte Carlo test (250 runs with randomized data). The coefficients of determination (r 2 ) between the original plot distances and distances in the final ordination solution were calculated to assess the community variability represented by the NMS axes [52]. The analysis was run using PC-Ord Software, version 7.03 [53]. Significant functional variables, CWM, were overlaid on the NMS ordination as correlation vectors [52] to confirm species classification. This was done looking at the distance between the position of the CWM of the functional groups and each species, taking also into consideration the position of the environmental variables.
Buffer distance for each land cover was selected individually, as its influence may vary with land cover type. Spearman correlations (to account for possible nonlinearity in the relationships) were performed between NMS scores, the CWM of each functional group, and the environmental variables. Correlations were made in R software [51] and considered significant for p < 0.05. The buffer distance selected for each land cover was the one presenting the higher number of significant correlations with the functional variables [54]. The CWM values of each functional group were then interpolated for the entire city of Lisbon using ordinary kriging. All spatial analysis was done in ArcMap (ESRI, version 10.5).

Results
In total, we identified 38 species (Table S2). Of these, two species had an original poleotolerance classification (Table 1) of "natural", comprising 5.26% of the total, and six species had a poleotolerance classification of "seminatural" (15.79%). The majority of the species (17) were classified as "low disturbance" (44.73%), and 13 species were classified as "high disturbance", representing the remaining 34.21% species.
The NMS ordination suggested a three-dimensional solution (minimum stress = 10.8%; P = 0.004), as increasing the number of dimensions did not significantly lowered the stress. The first axis accounted for 45.7% of the variance. The second axis explained 32.5% and the third explained 11.0% of the total variance (89.3%). As the third axis explains only a small percentage of the variance, the analysis focuses on the first two axes.

Confirming Species Classification
The first step was to overlay lichen species ordination scores in the NMS solution, colored by functional group, to understand if the poleotolerance ranking assigned to each species was correct ( Figure 2A).
The analysis of the ordination revealed five clearly misclassified species, i.e., located in the opposite side of the gradient where the CWM of that functional group is higher. Enterographa crassa (DC.) Fée was classified as natural and was reclassified to low disturbance. Lecanographa amylacea (Pers.) Egea and Torrente from natural to seminatural. Lepra amara (Ach.) Hafellner from high disturbance to low disturbance. Finally, Parmotrema hypoleucinum (J. Steiner) Hale and Ramalina lacera (With.) J.R. Laundon from seminatural to low disturbance. The CWM was recalculated using this new classification ( Figure 2B).

The Poleotolerance Response Trait as an Indicator of Urban Disturbance
Ordination scores and the CWM of all functional groups were correlated with the environmental variables to determine the size of the buffer to use for each land cover (Table S1). As previously stated, from the 25 land cover types, five were immediately discarded. Remaining land cover types were grouped based on their buffer distances and correlation trends. From there, four groups were created. The first group was the HDLC (high disturbance land covers) comprising five land covers (all at a 500 m buffer) associated with high human disturbance (discontinuous dense urban fabric; industrial, commercial, public, military, and private units; port areas; other roads and associated land; water). The remaining three groups comprise land covers associated with low human disturbance, the low disturbance land covers (LDLC), with three different buffer distances: LDLCa, LDLCb, and LDLCc, with 1000, 1500, and 2000 meter buffers, respectively. LDLCa comprises the arable land, herbaceous vegetation associations, and sports and leisure facilities land covers. LDLCb comprises the fast transit roads and associated land, land without current use, pastures, and permanent crops. LDLCc comprises the discontinuous very low-density urban fabric, forests, and isolated structures. The analysis of the ordination revealed five clearly misclassified species, i.e., located in the opposite side of the gradient where the CWM of that functional group is higher. Enterographa crassa (DC.) Fée was classified as natural and was reclassified to low disturbance. Lecanographa amylacea (Pers.) Egea and Torrente from natural to seminatural. Lepra amara (Ach.) Hafellner from high disturbance to low disturbance. Finally, Parmotrema hypoleucinum (J. Steiner) Hale and Ramalina lacera (With.) J.R. Laundon from seminatural to low disturbance. The CWM was recalculated using this new classification ( Figure 2B).

The Poleotolerance Response Trait as an Indicator of Urban Disturbance
Ordination scores and the CWM of all functional groups were correlated with the environmental variables to determine the size of the buffer to use for each land cover (Table S1). As previously stated, from the 25 land cover types, five were immediately discarded. Remaining land cover types were grouped based on their buffer distances and correlation trends. From there, four groups were created. The first group was the HDLC (high disturbance land covers) comprising five land covers (all at a 500 m buffer) associated with high human disturbance (discontinuous dense urban fabric; industrial, commercial, public, military, and private units; port areas; other roads and associated land; water). The remaining three groups comprise land covers associated with low human disturbance, the low disturbance land covers (LDLC), with three different buffer distances: LDLCa, LDLCb , and LDLCc, with 1000, 1500, and 2000 meter buffers, respectively. LDLCa comprises the arable land, herbaceous vegetation associations, and sports and leisure facilities land covers. LDLCb comprises the fast transit roads and associated land, land without current use, pastures, and permanent crops. LDLCc comprises the discontinuous very low-density urban fabric, forests, and isolated structures. Correlations between the CWM and the axes scores and the best buffers are shown in Table 2. The environmental variables that were significantly correlated with ordination axes scores were overlaid in the solution as correlational vectors ( Figure 2B). Correlations between site ordination scores and environmental variables suggest that the first axis (Table 2) is reflecting a gradient of anthropic disturbance. On the most disturbed side of the gradient, we have land covers associated with high anthropic presence and disturbance: Continuous urban fabric (200 m buffer) and land covers from the HDLC group. On the less disturbed side, we have land covers associated with seminatural areas comprising the LDLCa, LDLCb and LDLCc land-cover groups. The remaining land covers are not significantly correlated with the first axis. In addition, both NDVI values (NDVIgs and NDVIb100) and the area of the green spaces were also associated with the less disturbed side of the gradient. The second axis also revealed a gradient of disturbance, though less evident than the first: One side is associated with the discontinuous medium density urban fabric (1000 m buffer) and railways and associated land (2000 m buffer) land covers (both associated to higher disturbance levels); the other one with green urban areas (200 m buffer) land cover (associated to lower disturbance levels), the area of the green spaces, and the NDVIb100. The remaining land covers were not significantly correlated with this axis. For more detailed data, see Table S3). Table 2. Spearman correlation coefficients between lichen community ordination axes (NMS1, NMS2, and NMS3), lichens CWM (high disturbance, low disturbance, and seminatural) and the environmental variables. Significance level is indicated by the asterisk: * = p < 0.05; ** = p < 0.01; *** = p < 0.001. Environmental variables: Area is the area of each green space; NDVIgs is the NDVI value for the green space centroid; NDVIb100 is the NDVI value for a 100 m buffer around each green space centroid; CUF-Continuous urban fabric (S.L.: >80%); DLDUF-Discontinuous low density urban fabric (S.L.: 10-30%); DMDUF-Discontinuous medium density urban fabric (S.L.: 30-50%); GUA-Green urban areas; RAL-Railways and associated land; HDLC-High disturbed land covers; LDCLa-Low disturbed land covers (1000 m buffer); LDCLb-Low disturbed land covers (1500 m buffer); LDCLc-Low disturbed land covers (2000 m buffer). To understand if the poleotolerance response trait was, in fact, mediating lichen species' response to this disturbance gradient, the CWM of poleotolerance functional groups were also overlaid in the solution as vectors (Figure 3). The high disturbance functional group was associated with the most disturbed side of the gradient. The significant correlations between the high disturbance CWM values with several environmental variables corroborate this result ( Table 2). The low disturbance functional group was associated with the opposite side of the disturbance gradient ( Figure 2B). Again, the low disturbance result was corroborated by the correlations with the environmental variables (Table 2). Lastly, the seminatural functional group was also associated to the low disturbance side of the gradient, though to a lesser extent than the two other functional groups ( Table 2). This functional group was positioned in the middle of the first axis disturbance gradient but clearly separated by the other two functional groups in the second axis, being positioned on the higher disturbance side of the second axis gradient. However, this functional group was only significantly correlated with the discontinuous medium density urban fabric (1000 m buffer) and the LDLCb land cover group. Correlations between the CWM and the axes scores and the best buffers are shown in Table 2. The CWMs values of each functional group in the 41 sampling sites were interpolated for the entire city of Lisbon to have a spatial distribution of anthropic disturbance levels indicated by the poleotolerance trait. In general, spatial patterns of high disturbance and low disturbance functional groups indicate that the most disturbed areas in Lisbon (orange and red areas in Figure 4A; blue and green in Figure 4B) are located closer to the airport and the city center, while low anthropic areas can be found close to Monsanto and further away from the city center ( Figure 4A,B). The spatial pattern of the seminatural functional group shows no clear pattern in Lisbon, being diffuse throughout the south and southeast part of the city ( Figure 4C).   Table S4.

Discussion
This study presents a novelty regarding the application of the poleotolerance response trait as an ecological indicator of human disturbance in urban environments. This trait was able to differentiate highly disturbed green spaces from those with lower levels of anthropic disturbance. The two functional groups with higher tolerance to greater disturbed areas showed to be efficient in this task. On the other hand, the seminatural functional group also indicated low levels of disturbance but failed to indicate the green spaces more similar to seminatural habitats.
We were able to identify five misclassified species in the poleotolerance response trait (Enterographa crassa, Lecanographa amylacea, Lepra amara, Parmotrema hypoleucinum, and Ramalina lacera). These results confirm that lichens distribution is a result of multiple interacting biotic and abiotic factors that may result in different ecological behaviors regarding anthropic disturbance. The fact that this trait classification was elaborated based only on studies conducted in Italy may potentially increase the probability of misclassifying species from regions with contrasting biotic and abiotic factors (climate, for example). Such a situation has already been reported by Pinho et al. [55] for the trait classification regarding epiphytic lichens tolerance to nitrogen pollution. Other authors also found contrasting classifications in different regions from the United States of America, something they attributed to climatic differences [56].
According to Nimis [36], the poleotolerance response trait classification should reflect lichens' ability to occur in areas with different levels of human disturbance. The disturbance gradient along Lisbon green spaces has been recently attributed to a strong background pollution gradient [14]. Our results showed that the functional groups with high and low tolerance to disturbance were able to successfully identify such a disturbance gradient. We saw that green spaces with bigger areas, higher NDVI values, and surrounded by low disturbance land covers (LDLCa/b/c) were those with lower anthropic disturbance (high CWM of low disturbance functional group and low CWM of the high  Table S4.

Discussion
This study presents a novelty regarding the application of the poleotolerance response trait as an ecological indicator of human disturbance in urban environments. This trait was able to differentiate highly disturbed green spaces from those with lower levels of anthropic disturbance. The two functional groups with higher tolerance to greater disturbed areas showed to be efficient in this task. On the other hand, the seminatural functional group also indicated low levels of disturbance but failed to indicate the green spaces more similar to seminatural habitats.
We were able to identify five misclassified species in the poleotolerance response trait (Enterographa crassa, Lecanographa amylacea, Lepra amara, Parmotrema hypoleucinum, and Ramalina lacera). These results confirm that lichens distribution is a result of multiple interacting biotic and abiotic factors that may result in different ecological behaviors regarding anthropic disturbance. The fact that this trait classification was elaborated based only on studies conducted in Italy may potentially increase the probability of misclassifying species from regions with contrasting biotic and abiotic factors (climate, for example). Such a situation has already been reported by Pinho et al. [55] for the trait classification regarding epiphytic lichens tolerance to nitrogen pollution. Other authors also found contrasting classifications in different regions from the United States of America, something they attributed to climatic differences [56].
According to Nimis [36], the poleotolerance response trait classification should reflect lichens' ability to occur in areas with different levels of human disturbance. The disturbance gradient along Lisbon green spaces has been recently attributed to a strong background pollution gradient [14]. Our results showed that the functional groups with high and low tolerance to disturbance were able to successfully identify such a disturbance gradient. We saw that green spaces with bigger areas, higher NDVI values, and surrounded by low disturbance land covers (LDLC a/b/c ) were those with lower anthropic disturbance (high CWM of low disturbance functional group and low CWM of the high disturbance functional group). Only one study has previously used this trait in an urban context [23]. The authors saw that in Rome, low and high disturbance functional groups accounted for 51% of total species (38% and 13%, respectively), but no discussion of the relationship of its patterns and urban disturbance was made. In Lisbon, we found that 87% of the lichen species were classified in these two groups (55% and 32%, respectively), a number significantly higher than that found for Rome. Our results showed that spatially, the higher anthropic disturbance areas are located close to the airport and the city center ( Figure 4A,B), coinciding with those previously highlighted as with lower air quality [14]. Conversely, the spatial distribution of these functional groups showed that low disturbance green spaces are located in the area of Lisbon with the largest green space, Monsanto ( Figure 4A,B), also matching the area previously pointed as with best air quality [14]. Thus, our results show that high and low disturbance functional groups have not only the potential to be applied in urban environments, but they can also be a useful tool to spatially signal the disturbance areas in a city, which can be used as a tool for better management practices in cities.
Similarly to the low disturbance functional group, the seminatural functional group also indicated low levels of disturbance, although without a clear spatial pattern. This could be due to the fact that only a very small proportion of the sampled lichen species were representative of this functional group, and it was thus difficult to produce robust results with it. Only two species belonged initially to the natural functional group but were reclassified and merged with the seminatural. As a result, after the reclassification, there were five species classified as seminatural out of the 38 total (13.16%). The low representativeness of natural and seminatural species may be a clue pointing to the lack of truly natural habitats. We expected that the seminatural functional group would be significantly positively correlated with the area of the green spaces, as bigger areas would work as a stopper for a surrounding disturbance. However, that did not happen. We sampled five green spaces with areas equal or higher than 10 ha, expecting these to represent seminatural areas (Monsanto-1000 ha; Quinta das Conchas-23 ha; Parque Bensaúde-22 ha; Parque José Gomes Ferreira-10 ha, and Tapada das Necessidades-10 ha). Monsanto is by far the largest green space in Lisbon, with dense vegetation throughout the whole area. Quinta das Conchas and Parque José Gomes Ferreira are both close to the airport and to one of the busiest roads in Lisbon and both have sparse vegetation. Parque Bensaúde also has sparse vegetation, and it is located right next to another high traffic road (Avenida Lusíada). Lastly, Tapada das Necessidades is composed of dense vegetation and sits close to the port area and the river Tagus. Thus, though they are big enough to buffer disturbance effects, the fact that they have low-density vegetation and/or they are located next to high disturbance sites (like the airport or major roads) may somewhat prevent them from being true seminatural areas. This translates into a lack of robustness of results regarding this functional group. We should however note that Reference [23] concluded that the high percentage (~50%) of seminatural and natural poleotolerance classified epiphytic lichens, in Rome, was due to a successful city environmental planning. Thus, it is possible that other cities may have the necessary environmental characteristics to house this type of lichens species. The fact that, in our work, natural and seminatural functional groups had a low representation may be seen as a limitation for the application of this trait in highly disturbed cities. In fact, in Lisbon, it is very difficult, if not impossible, to identify and sample green spaces that fall into the description provided by Nimis et al. [36], characterizing the habitat for these two functional groups (Table 1). Thus, we recommend that future works should take this factor into consideration when preparing the sample design, so that in cities like Lisbon, with high disturbance levels, sampling sites are characterized as natural or seminatural areas, even if located in the outskirts of the urban perimeter are included.
In this work, we intentionally focused on a single trait, poleotolerance, with the purpose of understanding if this trait could be useful as an indicator of anthropic disturbances in an urban context. Our results showed a shift in lichen communities from low tolerance species to high tolerance species as we moved to more disturbed areas. A recent work [14] using the same methodology in the same study area explored three other, different traits, growth form, eutrophication tolerance, and humidity requirements, and did not observe any shifts in lichen communities mediated by these three traits. This is an indication that in urban environments, such as this encountered in Lisbon, the use of the poleotolerance trait may be more advantageous than using other, more commonly used traits, giving a clearer indication about the disturbance gradient. Humidity requirements, growth form, and eutrophication tolerance are frequently used to indicate disturbance related to the urban heat island effect (microclimatic variations) and nitrogen pollution. The fact that poleotolerance trait performed better here is the advantage of using it when we are not able to clearly disentangle those two environmental drivers, highlighting its potential as an integrated indicator of overall anthropic disturbance in urban settlements. The presence of pruina or secondary metabolites also showed promising results in South American cities as indicators of the effects of disturbances associated to urban environments [28,57]. In our work, we did not test these traits. The poleotolerance is, for now, restricted to Southern European cities that share some of the lichen flora of Italy. Pruina or secondary metabolites have the advantage of being independent of this classification, having therefore more potential to be applied outside Europe, so future works should try to explore it in Europe, to understand if they also have the potential to be used here. We also suggest, for future works focusing on this trait, that, if possible, datasets of population densities and distributions should be used, as these types of data can be a very suitable proxy for human disturbance levels inside the studied areas.

Conclusions
We conclude that the poleotolerance response trait can be used to distinguish areas with different anthropic disturbance levels in regions of Europe apart from Italy. Both the high disturbance and the low disturbance functional groups enabled us to identify a gradient of anthropic disturbance. However, the seminatural functional group did not, and that is probably related to the low number of species belonging to this functional group, possibly due to the lack of truly natural green spaces in our study area. Further studies should try to fully elucidate if this is a trait limitation or a sampling design issue. Because this is an integrative trait, it has the potential to be used even when we do not know the underlying environmental factors driving changes, signaling disturbed areas that should deserve further attention in future studies to individually assess the drivers of change, a very important tool under an urban management perspective.
To sum up, even with the limitations that may arise, mostly due to the fact that species classification may not be present in the Italian database, our work showed that the poleotolerance response trait has potential to be used in urban environments in other Southern European cities. Future studies should consider enlarging the trait database for Europe, allowing its use for wider geographical regions.

Conflicts of Interest:
The authors declare no conflict of interest.