Are Urban Communities in Successional Stasis? A Case Study on Epiphytic Lichen Communities

Urban areas may contain a wide range of potential habitats and environmental gradients and, given the many benefits to human health and well-being, there is a growing interest in maximizing their biodiversity potential. However, the ecological patterns and processes in urban areas are poorly understood. Using a widely applicable ecological survey method, we sampled epiphytic lichen communities, important bioindicators of atmospheric pollution, on host Quercus trees in urban parks of London, UK, to test if common patterns relating to lichen diversity are mirrored in urban green spaces. We found lichen diversity to be dependent on host species identity, and negatively related to local tree crowding. In addition, we found a strong negative effect of tree size on lichen diversity, leaving large trees as unexploited niches. A novel network analysis revealed the presence of only pioneer communities, showing the lichen communities are being held in successional stasis, likely due to the heritage effects of SO2 emissions and current nitrogen pollution and particulate emissions. Our study highlights that jointly assessing species richness, community structure and the successional stage can be key to understanding diversity patterns in urban ecosystems. Subsequently, this may help best determine the optimum conditions that will facilitate biodiversity increase within cities.


Introduction
Cities contain a diverse range of habitats and environmental gradients within relatively small areas [1]. They also present unique combinations of climatic, hydrological and biogeochemical conditions [1,2]. These factors, added to their proximity to centres of research and academia, make cities a promising study system for testing and developing ecological theories. Despite this, most ecology studies are still focused on non-urban areas. This means that many of the processes affecting urban biodiversity are poorly understood and simply applying results from non-urban systems may not always be appropriate [3].
Some ecological patterns appear to be the same in urban and rural ecosystems. Metapopulation theory, for example, has been effectively applied to both urban and rural systems and shows that factors such as patch connectivity [4], patch density [5] and patch area [6] affect urban metapopulations in the same way as non-urban ones. A test of the species area relationship (SAR) has also shown comparable results for urban and rural bird communities [7]. The same study also revealed similar latitudinal diversity gradients between urban and rural communities [7]. Seasonal shifts in agricultural fertilisers, remains high and is what now appears to be affecting current lichen assemblages in London and other cities most strongly [31,35]. London communities in particular are still low in biodiversity and are dominated by pollution-tolerant species such as Xanthoria parietina (L.) Th. Fr. 1860 and Physcia adscendens H. Olivier 1882 [36,37].
The effects of pollution on lichen diversity interact with many other ecological factors and understanding this interaction is key to elucidating patterns in urban lichen communities. For epiphytic lichen communities, one of the most reported factors in predicting diversity at the level of the individual tree is the trunk diameter, often measured as a diameter at breast height (DBH). DBH, which correlates with a number of other factors [38], and it affects lichen diversity partly because larger trees are also older and therefore may have had more time for lichen colonization [39]. However, from a pure ecological point of view a larger DBH implies more surface area for epiphytic lichens and so we can expect higher diversity from the classical species area relationship. Indeed, the predominant relationship of DBH with lichen (alpha) diversity has been positive across a range of study areas and host tree species [38][39][40][41][42][43][44][45][46][47]. However, examples with a unimodal [44,48,49], negative [50,51] and even no relationship at all also exist [28,[52][53][54][55], suggesting there may be no universal relationship between lichen diversity and DBH.
Functional categorisation can be a powerful tool in understanding responses to environmental factors. The use of functional traits, such as reproductive strategy, growth form and eutrophication tolerance, is a rapidly growing field of research in lichen ecology [25,26,56]. Individual lichen species can show strong preferences for bark chemistry and are often functionally categorised according to their substrate pH preference [31,57,58] and nitrogen tolerance [59,60]. Despite its widespread use in the literature, Giordani and Malaspina (2016) showed that under eutrophic conditions, bark pH alone did not significantly affect lichen community abundance [28]. This result therefore leads to the possibility that the effect of bark pH may be site-specific or dependent on the type of nitrogen source. Furthermore, bark pH also co-correlates with cation concentration and nitrogen pollution levels, making it difficult to separate the effects of these variables [14,35].
Bark structure and chemistry often differ between tree species [61], and consequently, lichens will preferentially associate with certain host tree species. The diversity of epiphyte communities can therefore be explained and predicted by the diversity of hosts [62,63]. It has actually been suggested that the tree species is of greater importance than bark pH alone in explaining lichen colonization, as a changing pollution climate concomitant with a general increase in bark pH seems to have made many lichen species less sensitive to substrate pH [64] as shown above in Giordani and Malaspina (2016).
The use of functional traits has also uncovered important successional trends that are not apparent when only considering species richness. Functional trait analyses along chronosequences of tree age reveal significant species turnover as trees age/grow, with a general shift from predominantly crust-like (crustose), sexually reproducing communities to leafier (foliose), asexual ones [49,[65][66][67]. This community succession can be explained by a combination of allogenic factors, relating to changes in microhabitat conditions as the tree gets older/larger, and autogenic factors such as out-competition of smaller early-stage pioneers by larger, late-stage species [14]. Different successional stages differ in their response to environmental factors. Therefore, investigating community composition and structure as well as species richness may help to explain certain patterns in lichen ecology.
Individual host trees as well as woodland fragments can be seen as potential habitat patches for lichen species, making epiphytic lichens an ideal study system to test metapopulation and metacommunity theory. Connectivity between individual trees and woodland patches could therefore be very important in determining colonisation of trees, which in turn is also dependent on the dispersal abilities of lichen species. Lichen spores and vegetative propagules are very small so we might expect them to be well dispersed by wind, yet various studies have uncovered evidence for dispersal limitation [14,38,[68][69][70][71], although most have focused on one or a few lichen species. For example, Buckley (2011) found that whilst tree DBH explained the most variation, lichen diversity was negatively related to distance to a source population of woodland in New Zealand mountain beech trees, indicating lichens are limited in their ability to colonize new substrate, even over distances of less than 1 km. Furthermore, Buckley also found a weak but significant negative effect of local crowding (tree density within 5 m of the host), and this could be caused by crowding leading to shading and a change in the conditions required for some lichen species to colonise, especially early succession species. Hence, epiphytic diversity can be affected by processes occurring at a number of spatial scales, and the tree host spatial ecological context could be important for determining biodiversity patterns.
The complex interaction between the abiotic and biotic factors discussed above makes it difficult to predict which have a significant effect on the diversity of urban epiphytic lichen communities. For this reason, we chose to investigate these factors simultaneously and model how their interaction may affect the diversity of epiphytic lichen communities in London. Specifically, we investigated how tree host species, its geographic location, spatial isolation and size/age, can be used to predict epiphytic Diversity 2020, 12, 330 4 of 18 lichen diversity in urban green spaces in London. We also looked at how functional trait diversity responds to these factors and whether determining the successional stage of the communities can help to understand the observed trends. We hypothesized that: 1.
The history of change in pollution levels in London would mean that older and therefore larger trees may have spanned both high and low atmospheric pollution and that this could alter the tree DBH-lichen diversity relationship.

2.
Lichen communities show a shift in functional trait diversity with sexual reproduction and basidophilous substrate pH preference becoming less abundant as tree size increases.

Survey Sites
We sampled 124 Quercus trees in Richmond Park, Regent's Park and Victoria Park in Greater London, UK ( Figure 1). These parks were chosen to get an approximate west to east transect covering a pollution gradient from lower in the west to higher in the east [72]. All three contained at least 40 Quercus host trees of at least two of the species included in this study (Quercus robur, Q. cerris, Q. rubra and Q. palustris), enabling us to test for interspecific variation of lichen communities within the genus Quercus. We selected trees of various sizes to cover a range of ages and successional stages of epiphytic communities. We georeferenced each tree using an iPhone 6s GPS, which has an approximate mean distance error of 4.19 m [73]. Individual tree locations can be found in Figure S1. diversity responds to these factors and whether determining the successional stage of the communities can help to understand the observed trends. We hypothesized that: 1. The history of change in pollution levels in London would mean that older and therefore larger trees may have spanned both high and low atmospheric pollution and that this could alter the tree DBH-lichen diversity relationship. 2. Lichen communities show a shift in functional trait diversity with sexual reproduction and basidophilous substrate pH preference becoming less abundant as tree size increases.

Survey Sites
We sampled 124 Quercus trees in Richmond Park, Regent's Park and Victoria Park in Greater London, UK ( Figure 1). These parks were chosen to get an approximate west to east transect covering a pollution gradient from lower in the west to higher in the east [72]. All three contained at least 40 Quercus host trees of at least two of the species included in this study (Quercus robur, Q. cerris, Q. rubra and Q. palustris), enabling us to test for interspecific variation of lichen communities within the genus Quercus. We selected trees of various sizes to cover a range of ages and successional stages of epiphytic communities. We georeferenced each tree using an iPhone 6s GPS, which has an approximate mean distance error of 4.19 m [73]. Individual tree locations can be found in Figure S1.

Host Species
We selected Quercus as the only host genus to reduce the variation in genus-specific traits (e.g., bark properties, nutrient content and moisture retention), and because being a widespread taxon it would facilitate comparison to previous studies and guarantee adequate sample sizes. The bark of Quercus in the UK tends to be slightly acidic, which provides a buffering effect with NH3 and/or NOx

Host Species
We selected Quercus as the only host genus to reduce the variation in genus-specific traits (e.g., bark properties, nutrient content and moisture retention), and because being a widespread taxon it would facilitate comparison to previous studies and guarantee adequate sample sizes. The bark of Quercus in the UK tends to be slightly acidic, which provides a buffering effect with NH 3 and/or NO x pollution, limiting nitrophytic lichens to truly eutrophic sites. High abundances of nitrophytic lichens on oak will hence indicate high concentrations of environmental nitrates [31] and these lichens can therefore be considered reliable bioindicators of urban pollution.
The species Quercus robur, Q. cerris, Q. rubra and Q. palustris were the most abundant in the three parks and all trees met the standardised EU lichen recording criteria: free-standing; over 70 cm in trunk girth; less than 10 • inclination from vertical and with no clear signs of disease/stress (for example burls or high Ganoderma coverage) [74]. We avoided trees near highly trafficked roads, where possible, to reduce localised effects of pollution. Q. robur is the only native UK host tree species in this study.

Sampling Method
We concentrated on lichen communities on the lower boles of trees and did not sample upper boles or twigs. Lichen biodiversity was sampled using the standardised EU lichen recording method; sampling quadrats were attached on the north, south, east and west faces of the tree at 1.5 m above ground level and the abundance of all lichen species within the quadrats was recorded (see [74] for full details of protocol). We measured the tree circumference at 1.5 m to give an estimation of tree age and size, and tree isolation was quantified via a number of proxies: (1) the distance to the nearest tree, (2) distance to the nearest Quercus and (3) the number of trees within a 10 m radius. We also estimated bark rugosity by separating bark structure into five classes following the quantitative scale of Öckinger 2005 [69].
We identified lichen species predominantly in the field using a 10× hand lens and C and K chemical spot tests following Dobson (2011) and the online key 'A key to common lichens on trees in England' [75,76]. Field identifications were confirmed in the lab using a Leica MZ6 dissecting scope (Leica, Wetzlar, Germany) and a Zeiss Primo Star confocal microscope (ZEISS, Oberkochen, Germany) to observe and measure spores and microscopic characters.

Statistical Analysis
We produced rarefaction curves ( Figure S2) using the 'vegan' package v.2.4-5 implemented in R v.3.6.2 [77] to make sure sampling intensity was sufficient. A lichen diversity value (LDV) was quantified for each tree by calculating the mean abundance (number of 10 cm × 10 cm quadrats occupied) across all lichen species at each cardinal point individually and then summating the four means (N, S, E and W) to get a single LDV per tree [74].

Generalized Linear Modeling
We used a generalized linear model (GLM) to test the effects of the environmental variables measured on lichen diversity (LDV), using a gamma distribution with a log link function. Gamma was used as LDV is a continuous variable bound at zero. The explanatory environmental variables used were tree circumference (continuous), host tree species (categorical with four levels), park (categorical with three levels), bark rugosity (categorical with five levels), distance to nearest tree (continuous), distance to nearest Quercus (continuous) and number of trees within 10 m (count data). We used the following interaction terms: host tree species × circumference, host tree species × bark rugosity and circumference × bark rugosity. Interaction terms were only included for those variables that could be clearly linked with a biological explanation; for example, we may expect the relationship between lichen biodiversity and tree circumference to differ depending on tree species. The GLM was implemented in R using the base function 'glm' [77]. We added a value of one to the response variable LDV as the glm function does not accept zero-values with a gamma distribution.
We used a top-down model selection approach, beginning with the full model, using partial F-tests to assess the significance of variables and interactions terms, removing non-significant variables in a stepwise fashion and testing the model again until we reached our minimum best-fitting model. Factor levels for categorical explanatory variables showing significant effects on LDV were collapsed Diversity 2020, 12, 330 6 of 18 into a single level if t-tests showed no significant difference between those factor levels. An ANOVA comparing the collapsed factor-level model and the non-collapsed model tested whether collapsing factor levels significantly affected the model. Residual vs. fitted plots and Q-Q plots accompanied each model to test for homoscedasticity and normality of errors. Generalized variance inflation factors were calculated to test for collinearity of variables.
For functional trait analysis, we assigned functional groups to each species for two traits, reproductive strategy (sexual vs. asexual-sorediate) and substrate pH preference (acidophilous vs. basidophilous). These were assigned according to [56,60] (see Table S1 for functional groups information). Non-binary functional traits were not included to simplify analyses. We then calculated the community weighted mean (CWM) of each functional group on each tree by weighting the functional group by the abundance of species belonging to that group [78]. This was computed using the 'dbFD' function in the FD package v.1.0-12 [79] in R. GLMs were repeated as with the LDV but using the CWM of each functional group as the output variable. A binomial distribution was used as the CWMs are proportions between zero and one. For binary functional traits (acidophilous vs. basidophilous and sexual vs. asexual-sorediate), we only produced models for one of the groups (acidophilous and sexual), as the CWM of the two groups of the trait is complementary, adding up to one.

Species Co-Occurrence and Community Detection
Finally, we assessed community structure using paired species co-occurrence probabilities, which can be used to detect species interactions and/or shared habitat preferences. Co-occurrences assessed whether lichen species were significantly positively, negatively or randomly associated with each other. We converted the abundance data into presence-absence to conduct co-occurrence analyses. Co-occurrence probabilities were calculated using the probabilistic model of Veech (2013) via the 'cooccur' package v.1.3 in R [80,81] with the threshold argument set to 'TRUE' to remove any species pairs expected to occur at less than one site. We used a 2.5% significance level for the p-values as the test is two-tailed (species can be positively or negatively co-occurring). We then used the effect sizes matrix produced by 'cooccur' for community detection using K-means clustering. The standardised effect size calculated by cooccur ranges from −1 to 1, where increasingly negative values indicate stronger negative co-occurrences and increasingly positive values indicate stronger positive co-occurrences. However, K-means clustering requires distance measures that must be positive. Therefore, effect sizes, E ij , for each lichen pair (i and j) were converted into a distance d ij using the relationship d ij = (1 − E ij ), so that positively co-occurring species pairs had low distance scores and negatively co-occurring species pairs had high distance scores. Species that showed no significant positive or negative associations from the cooccur analysis were removed from the distance matrix prior to running the clustering analysis. The data was partitioned into k = 2 clusters using 'partitioning around medoid (pam)' clustering using the pam function of the R package 'cluster' [82].
Raw data and computer code for all statistical analyses in R are available in Supplementary Materials.

Results
Overall, we identified 24 lichen species on Quercus trees across the three parks. The abundance for all species and grouped abundance of Physcia spp. is shown in Figure 2. Abundances of Physcia species were grouped as the thalli of most Physcia individuals (most likely P. adscendens and P. tenella) were immature and lacked the characters needed for species-level identification. All species have been previously recorded in London [36]. Table S1 gives the full list of species with assigned functional groups. The highest individual tree LDV was 3.708 and the lowest was zero. The mean LDV across all individual host trees was 1.239 with a standard deviation of 0.850.

Lichen Diversity Value (LDV)
The partial F-test of the full model showed bark rugosity, all of the proxies for tree isolation and all the interaction terms to be non-significant at a p < 0.05 significance level (Table S2). The nonsignificant terms were removed from the model in a stepwise fashion, removing the variable with the highest non-significant p-value and rerunning the model until only significant variables remained. This process showed that the variable Park was also non-significant. This resulted in a final model containing three significant explanatory variables: tree circumference, tree species and number of trees within 10 m. All three variables had a statistically significant effect on LDV independent of their order in the model. This was shown using partial F-tests, again at a p < 0.05 significance level (circumference: F1121 = 45.140, p = 6.84 × 10 −10 , species: F3117 = 13.570, p < 1.14 × 10 −7 , No. trees 10 m: F1121 = 8.044, p = 0.005).
There appeared to be no significant difference in LDV between some of the factor levels of tree species (see t-test p-values between Q. cerris (Intercept), Q. rubra and Q. palustris in Table S3). Q. cerris, Q. rubra and Q. palustris were therefore collapsed into a single factor level. An ANOVA showed no significant difference in the explanatory power between the collapsed model and the non-collapsed model (F118,120 = 0.118, p = 0.889, see Table S4 for a full ANOVA table). The collapsed model therefore became our minimum adequate model with tree species having two factor levels Quercus robur and Quercus cerris/rubra/palustris.
Diagnostic plots showed no violation of the model assumptions (see Figure S3) and variance inflation factors (vifs) showed no significant collinearity between variables (all vifs were below critical value of 5; circumference = 1.082, species = 1.035, number trees within 10 m = 1.078). The final model had a Nagelkerke-R 2 of 0.546 and a Cox and Snell R 2 of 0.496. Model statistics are shown in Table 1.

Lichen Diversity Value (LDV)
The partial F-test of the full model showed bark rugosity, all of the proxies for tree isolation and all the interaction terms to be non-significant at a p < 0.05 significance level (Table S2). The non-significant terms were removed from the model in a stepwise fashion, removing the variable with the highest non-significant p-value and rerunning the model until only significant variables remained. This process showed that the variable Park was also non-significant. This resulted in a final model containing three significant explanatory variables: tree circumference, tree species and number of trees within 10 m. All three variables had a statistically significant effect on LDV independent of their order in the model. This was shown using partial F-tests, again at a p < 0.05 significance level (circumference: F 1121 = 45.140, p = 6.84 × 10 −10 , species: F 3117 = 13.570, p < 1.14 × 10 −7 , No. trees 10 m: F 1121 = 8.044, p = 0.005).
There appeared to be no significant difference in LDV between some of the factor levels of tree species (see t-test p-values between Q. cerris (Intercept), Q. rubra and Q. palustris in Table S3). Q. cerris, Q. rubra and Q. palustris were therefore collapsed into a single factor level. An ANOVA showed no significant difference in the explanatory power between the collapsed model and the non-collapsed model (F 118,120 = 0.118, p = 0.889, see Table S4 for a full ANOVA table). The collapsed model therefore became our minimum adequate model with tree species having two factor levels Quercus robur and Quercus cerris/rubra/palustris.
Diagnostic plots showed no violation of the model assumptions (see Figure S3) and variance inflation factors (vifs) showed no significant collinearity between variables (all vifs were below critical value of 5; circumference = 1.082, species = 1.035, number trees within 10 m = 1.078). The final model had a Nagelkerke-R 2 of 0.546 and a Cox and Snell R 2 of 0.496. Model statistics are shown in Table 1.  Figure 3 visually depicts the results of the generalized linear model and shows that LDV declined significantly with increasing host tree circumference (Figure 3a) and number of trees within 10 m (Figure 3b). The same relationship was observed for all host species. Figure 3 also shows the significant effect of species with Q. robur supporting higher lichen diversity than Q. cerris/rubra/palustris.   Figure 3 visually depicts the results of the generalized linear model and shows that LDV declined significantly with increasing host tree circumference (Figure 3a) and number of trees within 10 m (Figure 3b). The same relationship was observed for all host species. Figure 3 also shows the significant effect of species with Q. robur supporting higher lichen diversity than Q. cerris/rubra/palustris.

. Scatterplots showing the relationship between epiphytic lichen diversity (LDV) and host
Quercus circumference, number of trees within 10 metres and Quercus species. Quercus cerris, Q. rubra and Q. palustris were collapsed into one factor level (blue) due to a lack of significant difference between them. (a) Lichen diversity (LDV + 1) significantly decreased with increasing tree circumference; (b) LDV significantly decreased with increasing number of trees within 10 metres. Tree species significantly affects LDV with Q. robur (red) hosting higher lichen diversity than the collapsed group of Q. cerris, Q. rubra and Q. palustris (blue). No significant interaction is seen between tree circumference, number of trees within 10 m and tree species. Fitted regression lines and 95% confidence intervals (coloured boxes around regression lines) were predicted using the minimum adequate model shown in Table 1.

Functional Trait Analysis
We found no significant effects between environmental variables and the CWMs of the functional traits reproduction type and pH preference (Tables S5 and S6).

Species Co-Occurrence and Community Detection
Of the 276 possible species pairs, 123 were removed as the expected co-occurrence fell below the cooccur threshold of 1 due to one or both species being very rare. Across the remaining 153 species Quercus circumference, number of trees within 10 metres and Quercus species. Quercus cerris, Q. rubra and Q. palustris were collapsed into one factor level (blue) due to a lack of significant difference between them. (a) Lichen diversity (LDV + 1) significantly decreased with increasing tree circumference; (b) LDV significantly decreased with increasing number of trees within 10 metres. Tree species significantly affects LDV with Q. robur (red) hosting higher lichen diversity than the collapsed group of Q. cerris, Q. rubra and Q. palustris (blue). No significant interaction is seen between tree circumference, number of trees within 10 m and tree species. Fitted regression lines and 95% confidence intervals (coloured boxes around regression lines) were predicted using the minimum adequate model shown in Table 1.

Functional Trait Analysis
We found no significant effects between environmental variables and the CWMs of the functional traits reproduction type and pH preference (Tables S5 and S6).

Species Co-Occurrence and Community Detection
Of the 276 possible species pairs, 123 were removed as the expected co-occurrence fell below the cooccur threshold of 1 due to one or both species being very rare. Across the remaining 153 species pairs we found 10 significant negative co-occurrences, 34 positive co-occurrences and 109 random co-occurrences (p < 0.025). These results are summarised in the species co-occurrence matrix Diversity 2020, 12, 330 9 of 18 in Figure 4a, which trims out any species with no significant positive or negative co-occurrences. From the results, two separate communities were identified visually with positive co-occurrences appearing within the communities and negative co-occurrences between the communities (Figure 4b).
The community on the left of Figure 4b will from here on be informally referred to as Parmelion, and the community on the right as Xanthorion. Partitioning around the medoid analysis supported the presence of these two communities when k was set to two (Figure 4c). Cluster 1 in Figure 4c (pink) matched the Xanthorion community of Figure 4b and cluster 2 (green) matched the Parmelion community. Dimension one of the cluster plot explained 47.4% of the variance and dimension two a further 13.9%.
Diversity 2020, 12, x FOR PEER REVIEW 9 of 18 pairs we found 10 significant negative co-occurrences, 34 positive co-occurrences and 109 random cooccurrences (p < 0.025). These results are summarised in the species co-occurrence matrix in Figure  4a, which trims out any species with no significant positive or negative co-occurrences. From the results, two separate communities were identified visually with positive co-occurrences appearing within the communities and negative co-occurrences between the communities (Figure 4b). The community on the left of Figure 4b will from here on be informally referred to as Parmelion, and the community on the right as Xanthorion. Partitioning around the medoid analysis supported the presence of these two communities when k was set to two (Figure 4c). Cluster 1 in Figure 4c (pink) matched the Xanthorion community of Figure 4b and cluster 2 (green) matched the Parmelion community. Dimension one of the cluster plot explained 47.4% of the variance and dimension two a further 13.9%.

Discussion
Here we focussed on lichen epiphytic diversity in a large, historically well-recorded capital city that has undergone an improvement in atmospheric pollution levels over the last few decades. Overall, and in agreement with most other studies, we found that tree species and tree circumference at breast height (a proxy for tree size and age) explained the majority of the variation in host tree lichen diversity. We found a significant negative relationship between tree circumference and lichen diversity, meaning large, older trees are a lesser used niche by lichens. Local tree crowding, measured via a number of trees within 10 m, also explained a significant proportion of the variation in lichen diversity, with less crowded trees hosting a higher diversity. In contrast, neither distance to the nearest tree, locality or bark rugosity could explain any significant variation in tree host lichen diversity. Finally, our species co-occurrence analyses led to the detection of two distinct lichen communities (Xanthorion and Parmelion), both of which relate to early-successional communities, and as we discuss below, these could be important indicators of why London park-based lichen community diversity was negatively related to host tree circumference.

Lichen Diversity Is Significantly Affected by Host Tree Circumference, Tree Species and Local Crowding
The relationship between host tree size and epiphytic lichen diversity is a complex one that seems to vary depending on the system under study. Examples of both linear positive [38][39][40][41][42]83] and negative correlations [50,51] have been reported, as well as no correlation at all [28,[52][53][54][55]. In some cases, a unimodal relationship has been observed, with lichen diversity increasing up to a point and then subsequently plateauing or even decreasing with very large tree sizes [44,48,49]. Various factors correlate with tree size including tree age, bark chemistry and structure, and crown size, which may obscure a clear relationship between tree size and epiphytic lichen richness [14]. The relationship may also depend on the tree species and age/size range being studied, as well as the climatic, edaphic and atmospheric conditions of the site. It is therefore difficult to generalise this relationship and demonstrates the importance of site-specific studies when investigating patterns in lichen diversity.
In urban systems, it is important to consider how atmospheric pollution may interact with the effect of tree size. Though current atmospheric sulphur dioxide (SO 2 ) levels in London are reasonably low, larger older trees have lived through past high pollution level environments, and this is still being reflected in their bark chemistry [27,84,85]. Such an effect is expected to be especially pronounced in long-lived, non-exfoliating trees like oak [85]. In addition, given that nitrogen pollution is still high in London [86], larger oaks may have had more time to accumulate these pollutants in their bark. The legacy effect of pollution on large trees has already been discussed in relation to initial lichen colonisation [27,31,84,85,87,88], with studies showing that following a drop in atmospheric pollution levels, small trees are often the first to be colonised by lichens. Therefore, the bark on the trunks of small oak trees in London may currently provide a chemically more favourable substrate for lichen establishment and growth than that of mature oaks. Whether substrate suitability in this system is due to pH, cations, heavy metals or nitrogen and other pollutants in the bark cannot be elucidated in this study and will require further testing. Besides pollutant effects, another difference between older and younger trees is that younger bark normally absorbs less water [89] and is more structurally stable/less likely to disintegrate [90]; both of these factors have been shown to benefit epiphytic lichen communities and may be contributing alongside the chemical properties to increase the suitability of young bark for lichen growth.
Light intensity is another factor that may contribute to the higher suitability of young trees for lichen growth. Larger trees tend to have denser canopies meaning that less light reaches the trunk [91,92]. Specialist, climax indicator lichen taxa, including Peltigera, Usnea and Lobaria species, are usually inhibited by high light intensities and are hence found predominantly in dense forests or on large trees [91]. Generalist, pollution-tolerant taxa such as species from the Xanthorion parietinae favour well-lit situations and are common in exposed areas [93]. As the species pool in this study comes largely from these more generalist communities, it follows that a higher diversity is found on the smaller trees where individuals receive more light. Light intensity may also explain the negative effect of local crowding we see in Figure 3b. Increased local crowding may reduce the amount of light reaching the tree boles, which could inhibit establishment and growth of the dominant London communities, which favour well-lit situations. It should be noted, however, that a recent study by Giordani and Malaspina (2016) using hierarchical partitioning of variance and generalized linear mixed models found light intensity had little to no effect on lichen functional diversity in a highly polluted system [28].
The significant effect of tree species on epiphytic lichen diversity has been observed previously [41,50,52,64,94]. The fact that we were able to see this in our study, using only four species of a single genus, shows how important tree species can be for understanding epiphytic lichen diversity. The differences in lichen diversity found amongst the Quercus species in our study could be due to various bark-related factors (pH, water retention and cation concentration) and further characterization of the microhabitats of each Quercus species is needed to understand why Q. robur hosts a higher diversity.
In contrast to tree host species, tree circumference and local crowding, we could not detect any effect of bark rugosity, tree isolation or geographic location (park) on lichen diversity. It is likely that the effects of rugosity are already included in the tree circumference and species factors, even though we did not find any evidence for its relevance even when these factors were not included in the analysis. The difficulty of separating the effects of bark-related variables such as rugosity, water-retention and bark thickness was highlighted by McDonald et al. (2017) who observed a similar pattern with little to no effect of bark fissuring on lichen community composition [57].
Whilst we did use a number of proxies for tree isolation, it may be the case that we did not estimate isolation at the best scale as we only considered very localized (within 10 m radius) isolation metrics surrounding the host tree. However, the parks are all relatively large and well populated with trees, and given lichen spores can travel long distances on the wind, it may be the case that dispersal limitation is not a significant factor, even in a large urban site. None-the-less future work should investigate the population genetic structure of the communities, both within and between parks as this might uncover more information about the level of lichen dispersal in an urban context. Recent research combining genetic information into a hierarchical Bayesian model has concluded that, at least for Lobaria pulmonaria, a species found in old growth forests, spore dispersal may occur over several hundred metres to kilometres [95].
The lack of an effect of a park can be explained by previous work from Larsen et al. (2007). In an attempt to study the lichen and bryophyte distribution on oak in Central London in relation to air pollution and bark acidity, the authors identified four distinct biogeographic zones for the different epiphyte communities based on lichen and bryophyte frequencies [31]. The parks included in our study all fall within the same zone (Group D2) according to Larsen et al. and should therefore not be expected to differ significantly in their diversity and species composition.
Finally, contrary to our hypothesis 2 in the introduction, we found no significant effects of host tree-related factors on the functional trait diversity measures. This was unexpected as previous studies have shown strong correlations between lichen functional traits and tree-related factors such host age/size; for example, lichen epiphyte communities on various trees species in Scotland show a clear shift from sexual crustose communities to asexual foliose ones as tree size/age increases [49,[65][66][67]. Species diversity in this study is considerably lower than in the previously mentioned studies, which in turn signifies lower functional trait diversity. This lack of variation in functional traits suggests that all the lichen species may be occupying a similar niche (discussed more below) and the functional groups that would benefit from tree-related factors like tree size are simply not present. Giordani and Malaspina (2016) observed a similar pattern in a highly polluted area of northwest Italy where the abundance of nitrophytic and oligotrophic functional groups showed little to no response to tree-related environmental predictors like tree size, bark pH and light. The authors showed that functional group abundance was much more strongly related to levels of atmospheric pollutants, mainly nitrogen compounds.

Species Co-Occurrence Reveals Two Distinct Early Successional Lichen Communities
Although the search for species co-occurrence patterns has a long history in ecology [96], it is only recently that they have been used to infer subcommunities, which may relate to key habitat types [97] or possibly successional stages. The two lichen communities that emerge from the co-occurrence analysis show clear negative associations between each other suggesting they are indeed separate communities, each requiring slightly different environmental conditions. Previous work on UK lichen communities supports the two communities presented here: Parmelion and Xanthorion. The Xanthorion community (right of Figure 4b, cluster 1 in Figure 4c) shares four of its seven species with the Physcietum ascendentis association within the Xanthorion parietinae alliance described by James, Hawksworth and Rose (1977) [93]. Physcietum ascendentis consists of species adapted to nutrient-enriched substrates and high light intensities [93]. The Parmelion community (left of Figure 4b, cluster 2 in Figure 4c) shares six of its eleven species with the Parmelietum revolutae association within the Parmelion perlatae alliance [93]. This alliance occurs mainly on well-lit bark of Fraxinus, Quercus and Larix trees and appears more pollution-sensitive than the Physcietum ascendentis [93]. The Xanthorion and Parmelion communities discussed here also overlap with the communities described by Ellis et al. (2015). The species within Xanthorion most frequently occur in Ellis et al.'s community A, a pioneer community found on smooth-barked young trees. In contrast, the species within Parmelion occur most frequently in their community D, another early successional community associated with mesotrophic habitats. Ellis et al.'s communities A and D have also been linked to James, Hawksworth and Rose's Xanthorion parietinae and Parmelion perlatae, respectively [98]. Our community analysis shows that despite vast improvements in the pollution level in London, the lichen biota is still dominated by pollution-tolerant communities.
Our co-occurrence analysis therefore only uncovers the generalist, early successional lichen communities. Results suggest the negative relationship observed between tree DBH and lichen diversity may, in part, be driven by a preference of these communities for the well-lit, exposed substrates of young trees and an absence of late-successional epiphytes that tend to occupy the niches of larger trees. The lack of late-successional lichen species, and their associated functional traits, is also likely to be an important factor behind our inability to detect a relationship between functional diversity and any of the variables we sampled.
Large trees provide important microhabitats for many lichen species and are important refuges for sensitive climax communities [99,100], yet we believe the current crop of older and larger trees in London could be carrying the legacy of previously high levels SO 2 pollution and have had a longer exposure time to nitrogen pollution, which still remains high. It is also possible that the trees are indeed suitable for colonisation by late-successional species, but that dispersal in from regions outside London has not yet occurred. However, we feel this is unlikely given recent evidence for long-distance dispersal [95].

Conclusions
Lichens make up a large proportion of the diversity of many ecosystems and form an integral part of trophic webs [14]. Understanding the factors that control their diversity is therefore vital in the conservation of all species. The ecology and conservation of epiphytic lichen communities is especially relevant in cities where green spaces rely heavily on trees to encourage wildlife. This study shows that assessing community structure and successional stage can be key to understanding diversity patterns in urban ecosystems. This will subsequently help best determine the optimum conditions that may facilitate biodiversity increase within cities. Though recent data suggest average occupancy of UK lichens has increased since 1970 [101], atmospheric pollution needs to be targeted in order for urban communities to recover to a preindustrial state. Due to the heritage effects of SO 2 emissions and current nitrogen pollution and particulates it may be many years before the bark of larger trees becomes hospitable for lichens and urban lichen communities may remain at a pioneer stage. This successional stasis limits the overall species and functional diversity, which may also lead to more vulnerable communities (a pattern that has already been observed in Italian lichen communities [102]). Recovery of London's epiphytic lichen communities is therefore dependent upon decreasing atmospheric pollutants followed by potentially subsequent dispersal events from neighbouring rural communities. Until this occurs, we expect trees with suitable sizes for more stable communities to remain devoid of lichens and their high potential microhabitats will remain empty. As it stands, these trees are unexploited niches.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/9/330/s1, Figure S1: Maps showing the locations of the trees sampled in this study, Figure S2: Rarefaction curves showing the accumulative number of species found when increasing the number of sites (trees) sampled, Figure S3: Diagnostic plots of minimum adequate model LDV~Circumference + Species + Number trees within 10 m, Table S1: List of identified species in this study including their assigned species number by the British Lichen Society and functional trait groups, Table S2: Results of partial F-test of generalized linear model: LDV + 1 = Circumference + Species + Park + Species × Circumference + Species × Rugosity + Circumference × Rugosity + Distance to nearest tree + Distance to nearest Quercus + Number trees within 10 m, Table S3: Summary table of generalized linear model: LDV + 1 = Circumference + Species + Number trees within 10 m, Table S4: ANOVA table comparing the non-collapsed (Model 1) and collapsed (Model 2) generalized linear model, Table S5: Results of χ 2 test of generalized linear model: CWM sexual reproduction = Circumference + Species + Park + Distance to nearest tree + Distance to nearest Quercus + Number trees within 10 m + Species × Rugosity + Circumference × Rugosity, Table S6: Results of χ 2 test of generalized linear model: CWM acidophilous pH preference = Circumference + Species + Park + Distance to nearest tree + Distance to nearest Quercus + Number trees within 10 m + Species × Rugosity + Circumference × Rugosity, Dataset 1: Measured host tree variables, lichen species abundances, and Lichen Diversity Values (LDV) for each sampled tree, Dataset 2: Functional trait groups for all lichen species identified in this study, Computer Code: Computer code for all statistical analyses in R.