Factors Influencing Epiphytic Lichen Species Distribution in a Managed Mediterranean Pinus nigra Arnold Forest

Lichens have important ecological functions in black pine forests, such as nitrogen fixation and nutrient cycling. Understanding lichen diversity could provide a better understanding of black pine ecosystems. The aim of this study was to identify the factors affecting the composition of lichen communities and their specific diversity in Mediterranean black pine forests. Research was conducted in 48 sampling plots. For the analysis, presence–absence and frequency data of lichen species were used. For stand level analysis, four community composition tables were created. We used bioclimate, topography, stand, and parent rock as variables. A total of 33 epiphytic lichen species were identified in the black pine forests from 282 sampled trees. Indicator lichen species were determined according to geographic region and stand age classes. Hypocenomyce scalaris was found to be an indicator species for old forests. Frequency data were more useful for revealing lichen species composition than presence–absence data. Of the topographic variables, elevation was the most prominent and had the highest explanation ratio for the composition of lichen species with a coefficient of correlation (R2) value of 0.49. Significantly positive (p < 0.001) relationships were found between epiphytic lichen richness and tree crown height, tree height, and bark pH. Our results revealed that to retain the trees in the stands rich in lichen species diversity is recommended in the managed forests.


Introduction
Mediterranean ecosystems have extremely high species diversity that can be managed within a conservation area [1].Black pine tree with its five subspecies, is distributed in European and Asian mountainous areas [2][3][4].The greatest distribution of black pine was recorded in Anatolia (Turkey), accountings for 2.5 million ha of its total worldwide distribution of 3.5 million ha [5].The most commonly distributed altitude is 300-1800 m above sea level (a.s.l.) in sub-Mediterranean and Mediterranean-montane regions [6].One of the optimum growing sites for black pine in Turkey is the Alaçam mountains [7] where rarely encountered double-layered stands of these trees can be found [5].The Alaçam mountains have been found to have substantially high plant biodiversity [8], which is attributed to their elevation.The Alaçam mountains have a more continental climate due to their geographical position compared to Mt. Olympus and Mt.Trapezitsa (Greece).Pinus nigra Arnold is a common tree species in Mediterranean ecosystems as well as in these mountains, and the trees are known for hosting a vast number of lichen species.
The distribution of lichen species and the factors affecting their distribution must be determined to understand the functions of lichens.Lichen diversity is influenced by both the host tree species (tree height, age, diameter, branch density, pH of tree bark, bark structure, and water holding capacity of bark), by the position of trees in the stand (canopy closure, species mixture, number of trees per ha), and the site quality (aspect, slope, and altitude) [29][30][31][32][33][34][35][36][37][38][39].
Topographical variables influence the water and nutrient budget of a given site [40].Many topographical variables have been included in models to explain the plant species distribution and plant society [40][41][42][43][44].To determine the effects on lichen species diversity, some limited topographical variables are used, like altitude, exposure, and slope.The northness, general curvature, catchment slope, and slope length have rarely been incorporated in studies seeking to understand the relationship between lichen species diversity and topographical variables [45].
Burgaz et al. [46] reported that altitude, slope, tree diameter, and tree height affect epiphytic lichen species distribution in Spain.Lichen species diversity changes due to tree species in Mediterranean forests [19,47,48].Lichen species distribution has been studied on the plane [49], oak [15], and beech [50] species in the Mediterranean.Belinchón et al. [15] found that lichen species vary according to the spatial position of the semi-deciduous forest trees species, i.e., whether they are located at the edge or the center of the stand.Besides the location of the trees, the distance of trees to the streams [51] and management strategies [17] affect the distribution of lichen species in the Mediterranean forests.
Several studies examined lichen species and their distribution in black pine forests [52][53][54].Pirintsos et al. [52] detected 23 lichen species at Mount Olympus and 63 species were recorded at Mount Trapezitsa, both in Greece [53].Guvenc et al. [54], recorded 20 species on black pine trees on Uluda g Mountain in Turkey.
In the present study, the distribution of lichen species was investigated at the stand level, and factors affecting the distribution were categorized into groups; climate, topography, stand, and bedrock types.We investigated the relationship of the richness of the lichen species with tree size attributes.Within the current study, we intended to clarify whether (i) bioclimatic, topographic, stand, and parent rock variables influence the lichen species composition; (ii) there is a particular lichen species that acts as an indicator for a geographic region or a life stage of black pine stands; and (iii) black pine tree size and bark properties affect lichen species richness.

Study Area
The study area included the Ulus, E grigöz, and Alaçam mountains, which together form the Alaçam Mountain range.The Alaçam mountain range is located partly in the cities of Balıkesir and Kütahya.The study site is located 29 • 15 30-28 • 15 00 E and 39 • 38 00-39 • 07 30 N (Figure 1).Based on the meteorological data gathered from the vicinity of study area, remarkable severe summer droughts occur in the region.The general climate type in the study area is Mediterranean, with the lowest and highest annual precipitation with 458.2 and 860.3 mm in Tavşanlı and Simav towns, respectively [55].
The most abundantly parent rocks, the original rock from which the soil was formed, are granit, tuff-anglomera, mélange, quartzite and dacite.[55].The most abundantly parent rocks, the original rock from which the soil was formed, are granit, tuff-anglomera, mélange, quartzite and dacite.The single-layer Anatolia black pine (Pinus nigra Arnold subsp.pallasiana (Lamb.)Holmboe) forests in the study area are scattered over an area of 91,744 ha [55].Some parts of the black pine forests are degraded, and a greater portion of the forests are productive forests managed by the age class method [56].Black pine forests regenerate through the large area shield method [57].The silvicultural maintenance period varies between 5 and 10 years, and the removed stem bark remains in the stand.

Species Data
Epiphytic lichen samples were collected from 95 long-term observation sample plots.All epiphytic lichens were recorded from August to November 2007 within sampled plots.At each plot, five single stems were randomly selected, and the lichen samples were collected from lower trunks (ranging from 0 to 2 m above the ground) of living black pine trees.
The collected epiphytic lichens were treated with chemical spot test and diagnosed with stereomicroscope [58][59][60].The distribution of species into morphological structures was as follows: 7 foliose, 15 fruticose, 10 crustose, and 1 squamulose.To eliminate the spatial autocorrelation, the sampling plots within 1 km or less from each other were discarded randomly with spThin package [61] in R (R Foundation for Statistical Computing, Vienna, Australia).The autocorrelation threshold was set to 1 km in a study conducted in beech forests harboring a high lichen species diversity in Bayern, Germany [62].Since the resolution distance of bioclim variables was 1 km, and also regarding the previous study, we applied the autocorrelation threshold value of 1 km.The remaining 48 sampling plots (Dursunbey, n = 26; Sındırgı, n = 16; and Simav, n = 6) were subjected to stand level The single-layer Anatolia black pine (Pinus nigra Arnold subsp.pallasiana (Lamb.)Holmboe) forests in the study area are scattered over an area of 91,744 ha [55].Some parts of the black pine forests are degraded, and a greater portion of the forests are productive forests managed by the age class method [56].Black pine forests regenerate through the large area shield method [57].The silvicultural maintenance period varies between 5 and 10 years, and the removed stem bark remains in the stand.

Species Data
Epiphytic lichen samples were collected from 95 long-term observation sample plots.All epiphytic lichens were recorded from August to November 2007 within sampled plots.At each plot, five single stems were randomly selected, and the lichen samples were collected from lower trunks (ranging from 0 to 2 m above the ground) of living black pine trees.
The collected epiphytic lichens were treated with chemical spot test and diagnosed with stereomicroscope [58][59][60].The distribution of species into morphological structures was as follows: 7 foliose, 15 fruticose, 10 crustose, and 1 squamulose.To eliminate the spatial autocorrelation, the sampling plots within 1 km or less from each other were discarded randomly with spThin package [61] in R (R Foundation for Statistical Computing, Vienna, Australia).The autocorrelation threshold was set to 1 km in a study conducted in beech forests harboring a high lichen species diversity in Bayern, Germany [62].Since the resolution distance of bioclim variables was 1 km, and also regarding the previous study, we applied the autocorrelation threshold value of 1 km.The remaining 48 sampling plots (Dursunbey, n = 26; Sındırgı, n = 16; and Simav, n = 6) were subjected to stand level analyses in our study (Figure 1).Therefore, the species found in the 37 discarded sampling plots, such as Alectoria sarmentosa, Cladonia fimbriata, Cladonia ramulosa, and Parmelia sulcate, were excluded.In the 48 plots, lichen data were collected from a total of 282 trees.The lichen diversity was described using two univariate variables (species presence-absence and species frequency) and one multivariate (species community composition) variable.The presence of epiphytic lichen species in five trees in each sampling plot was converted into the frequency of lichens.A species community composition table was created using both presence-absence (LPA) and frequency (LFQ) data.The most abundant species covered in a portion of 30% were determined using the ordiselect function in the goeveg package [63].Those selected species included: Bryoria capillaris, Buellia griseovirens, Hypocenomyce scalaris, Hypogymnia farinacea, Hypogymnia physodes, Hypogymnia tubulosa, Ochrolechia turneri, Platismatia glauca, Pseudevernia furfuracea, and Usnea subfloridana.The dataset of the most abundant 30% of species selected from LPA was renamed GLPA and that from LFQ was renamed GLFQ.Thus, four community composition data tables were created and subjected to stand-based analyses.Lichen species richness, tree dimensions, and bark properties were used for tree level analyses in all sampling plots (n = 95).After the missing tree height and canopy thickness data had been discarded, the 373 remaining trees were used in analyses.For tree-level analyses, the number of epiphytic lichen species per tree was used as the response variable.

Environment Data
Environmental variables were grouped into four main topics: bioclimatic, stand structure, topography and parent rock data (Table 1 and Appendix A).The variables for environment were sourced from a 30 arc-second (0.00833 • , under 1 km) resolution dataset of 19 bioclimatic variables from WorldClim 1.4 [64] and the ENVIREM [65] dataset, in which an additional 16 climatic and 2 geomorphological variables were included.Besides the climatic data, we used the microscale topographic variables calculated from field analyses, which also identify the ecological processes.SAGA GIS terrain analysis functions were used, and a total of 40 variables (slope, aspect, and curvature) were obtained from ASTER DEM at a 30-m resolution [66].GRASS GIS was used to calculate the northness and eastness of the data [67].
The parent rock types, the only categorical variable in this study, were gathered from 1:25,000 scaled geological map provided by Mining Detection and Seeking General Directory (MTA).
In this study, forest structure variables affecting the lichen species richness and composition in Mediterranean forests were used.The ages of trees were determined by counting the year rings on the pencils cored from four to seven trees per stand.The canopy closure values were predicted via visual estimation, ranging between 0.1 and 1.0.The crown height refers to the height from ground level to the lowest level of the crown, and dry branch thickness is the distance between the lowest and highest dry branches on the black pine stem.Both crown height and dry branch thickness were measured using a VERTEX III Laser (Langsele, Sweden) instrument.The tree diameter at breast height (Dbh) were measured using calipers in two perpendicular dimensions.
The pH of black pine bark was measured from a 1:10 diluted deionized (DI) water solution of ground bark [30,68,69].Total nitrogen content of the barks was determined using the Kjheldahl method [70] with a Buchi Auto Kjheldahl Unit K-370 device.

Statistical Analyses
Permutational multivariate analysis of variance (PERMANOVA) tests based on the Bray-Curtis ecological distance similarity indices were used to test the statistical significance of the epiphytic lichen community composition variances between study regions, age class (young age, <60 years; mature age, 60-120 years; and old age, >120 years), and parent rock types (granite (g), tuff-aglomera (ta), mélange (m), quartzite (q), and dacite (d)) using the "Adonis" function in the vegan package of R [71].
Indicator species analysis [72,73] was used to explain the relationship between the epiphytic lichen species and study regions and the age class of stands using the multipatt-multi-level pattern analysis function of the indicspecies R package [74].Indicator species for one region and one age class or combination of region and age classes were calculated.Based on the indicator species analysis, the indicator species for one cluster solely, or combinations of clusters, were calculated.Only significant species at the 0.05 level with IndVal values greater than 0.5 were assessed further.
Ordination methods were used to investigate the relationships between epiphytic lichen species composition and environmental variables.To compare the explanatory power of various environmental variable sets, the variations of the four Y response tables operated with the varpart function in vegan were classified [71].Partial redundancy analysis (RDA) was used for community matrices, and a partial multiple regression analysis was performed for vector-independent variables.
For the final RDA model, individual explanatory variables were used (Table 1, Appendix A).The significant explanatory variables were selected using the ordistep function in vegan [71].Before variation partitioning and statistical selection, multicollinearity between the explanatory variables was checked by pairwise correlations [38].Strongly correlated variables (r > 0.75) were excluded from the selection.The effect of explanatory variables was tested using F-statistics via Monte Carlo simulation with 999 permutations.The accepted significance level was 0.05.
The 8 most frequently found lichen species were detected in 48 investigated sample plots [54].The variables were determined as being the most frequently encountered in the lichen species composition LPA dataset.Subsequently, species response curves of the most abundant lichen species were drawn using generalized linear models (GLM) [63].
The Pearson pairwise correlation method was applied for the tree-level analyses and significance values are indicated [75].A linear regression analysis was applied to determine the lichen species richness and single tree properties, like the diameter at breast height, the tree height, crown height, dry branch thickness, the bark pH, and the N content of bark.All statistical analyses were conducted under the framework of R, version 3.5.2(R Foundation for Statistical Computing, Vienna, Australia, 20 December 2018)-Eggshell Igloo [76].
The analysis of the variance of the lichen species composition, including three regions, five main parent rock types, and three age groups, was tested with a permutational multivariate analysis of variance using the Bray-Curtis distance matrix.The analysis showed that the community composition of the lichens was significantly different in terms of compositional similarity (Bray-Curtis, R 2 ADONIS = 0.15, p = 0.001) in different regions, but no significant difference was detected for either parent rocks (Bray-Curtis, R 2 ADONIS = 0.11, p = 0.142) nor for age groups (Bray-Curtis, R 2 ADONIS = 0.06, p = 0.204).The indicator species analysis of the LPA data only identified two species as significant indicator species in Simav: Bryoria implexa and Letharia vulpina.Platismatia glauca was identified as the only indicator species for the Dursunbey and Sındırgı regions (Table 2).The LFQ data indicator species analysis identified Usnea subfloridana as the only indicator species for the Simav region (Table 2).Hypocenomyce scalaris was determined as the indicator of old forests.No specific lichen species was assigned as a significant indicator for young or mature forests or for parent rock groups.

Stand-Level Lichen Species Composition
The climate, topographic, and stand structure variables explained the largest part of the species compositional variance (Figure 2).Frequency data sets explained the lichen species composition more efficiently compared to the presence-absence dataset (Figure 2).The best explanation variation partition score of 0.72 was produced by the GLFQ data; LFQ, LPA, and GLPA resulted in explanation levels of 0.56, 0.34, and 0.26, respectively.Within the GLFQ dataset, topographical variables had an explanation level of 0.49; bioclimatic, parent rock, and stand structure variables had explanation rates of 0.32, 0.30, and 0.15, respectively (Figure 2).The corrected explanation power values for four datasets were: GLPA R 2 = 0.31, GLFQ R 2 = 0.25, LPA R 2 = 0.25, and LFQ R 2 = 0.23 (Table 3).Six topographical variables, five stand variables, and three bioclimatic variables were incorporated into the analyses (Table 3).The altitude and mean diameter at breast height were included in all analyses.Bio15 and mCH were included in three analyses, whereas the remaining variables were included in one or two analyses (Table 3).
After forward selection, 10 explanatory variables for LPA, 7 explanatory variables for GLPA, 7 explanatory variables for LFQ, and 6 explanatory variables for GLPQ were detected as significant in the final RDA model (Table 3).Among the topographical variables, only the elevation (Elev), catchment slope (CtchSlp), downslope curvature (DwnCurv), hill index (HI), tangential curvature (TnCurv), and texture (Txtr) were explanatory for the lichen species composition (Table 3).Stand variables, such as the mean diameter breast height (mDbh), mean crown height (mCH), total basal area of black pine in ha −1 (hecTBA), maximum height (maxH), and standard deviation of crown height (sdCH) revealed significant relationships with the lichen species composition (Table 3).Only three of the bioclimatic variables-precipitation seasonality (Bio15), index of the degree of water deficit below water need (AridxThorn), and count of the number of months with a mean temperature greater than 10 • C (MonthTemp10) (Table 3)-were found to have a significant relationship with the lichen species composition.Parent rock was also found to have a significant relationship with the lichen species composition (Table 3).As a result of the RDA between the LPA data set and environmental variables, a strong and positive relationship with elevation on the first axis of RDA and a negative correlation with AridxThorn, MonthTemp10, and maxH were indicated.H. scalaris was shown to be correlated with AridxThorn and MonthTemp10.T. chlorophylla and H. physodes showed negative correlations with maxH.The second axis of RDA was negatively correlated with Txtr and TnCurv.U. subfloriana and L. vulpina were negatively correlated with Txtr.B. griseovirens was negatively correlated with hecTBA.The corrected explanation power values for four datasets were: GLPA R 2 = 0.31, GLFQ R 2 = 0.25, LPA R 2 = 0.25, and LFQ R 2 = 0.23 (Table 3).Six topographical variables, five stand variables, and three bioclimatic variables were incorporated into the analyses (Table 3).The altitude and mean diameter at breast height were included in all analyses.Bio15 and mCH were included in three analyses, whereas the remaining variables were included in one or two analyses (Table 3).
After forward selection, 10 explanatory variables for LPA, 7 explanatory variables for GLPA, 7 explanatory variables for LFQ, and 6 explanatory variables for GLPQ were detected as significant in the final RDA model (Table 3).Among the topographical variables, only the elevation (Elev), catchment slope (CtchSlp), downslope curvature (DwnCurv), hill index (HI), tangential curvature (TnCurv), and texture (Txtr) were explanatory for the lichen species composition (Table 3).Stand variables, such as the mean diameter breast height (mDbh), mean crown height (mCH), total basal area of black pine in ha -1 (hecTBA), maximum height (maxH), and standard deviation of crown height (sdCH) revealed significant relationships with the lichen species composition (Table 3).Only three of the bioclimatic variables-precipitation seasonality (Bio15), index of the degree of water deficit below water need (AridxThorn), and count of the number of months with a mean temperature greater than 10 °C (MonthTemp10) (Table 3)-were found to have a significant relationship with the lichen species composition.Parent rock was also found to have a significant relationship with the lichen species composition (Table 3).As a result of the RDA between the LPA data set and environmental variables, a strong and positive relationship with elevation on the first axis of RDA and a negative correlation with AridxThorn, MonthTemp10, and maxH were indicated.H. scalaris was shown to be correlated According to our results for the RDA between the GLPA data set and environment variables, the first axis of RDA is positively correlated with elevation.The second axis showed a negative correlation with mCH.We observed that H. scalaris and H. physodes species were affected negatively by mDbh together with the climate (Figure 3).
For the RDA performed between LFQ data and environment variables, the first axis of the RDA showed a positive correlation with elevation and a negative correlation with HI (Figure 3).
The first RDA axis indicated a positive correlation with mCH as a result of the RDA performed between GLFQ data and environmental variables.The second axis showed a positive correlation with Bio15.H. physodes was correlated with HI (Figure 3).
A high possibility of the presence of Pseudoverniea furfuracea, Hypogymnia farinacea, and Hypogymnia tubulosa was found regarding their abundance distribution (Figure 4A-F).We estimated that Plastimatia glauca was not present above 1450 m, and Hypogymnia physodes was not present above 1100 m a.s.l.(Figure 4A).The higher the average stand diameter at breast height (mDbh) was, the higher the probability of the presence of Bryoria capillaris and H. physodes; for Platismatia glauca, Ochrolecha turneri, and Buella griseovirens, the probability decreased (Figure 4B).The possibility of the presence of B. capillaris decreased and increased for P. glauca, as long as the increased precipitation seasonality (coefficient of variation) value of Bio15 increased (Figure 4C).Typically, the presence of B. capillaris and B. griseovirens decreased with increasing AridxThorn (Figure 4D).The analyses predicted that B. capillaris and B. griseovirens may disappear above the AridxThorn values of 72 and 68, respectively (Figure 4D).B. griseovirens and O. turneri were present above an altitude of 1100 m a.s.l.The probability of the presence of B. capillaris and H. physodes increased with increasing mCH (Figure 4E).O. turneri, P. glauca, B. capillaris, and B. griseovirens were present with a higher probability with a higher CtchSlp value (Figure 4F).

Tree-Level Lichen Species Richness
A Pearson pairwise correlation analysis was conducted between epiphytic lichen species richness and individual tree properties, including the Dbh, dry branch thickness, crown height, total nitrogen, and bark pH (Figure 5).Correlations between lichen species richness and Dbh and total nitrogen of bark were not significant (Figure 5).We detected positive correlations between the epiphytic lichen richness and crown height and tree height; and negative correlations with bark pH were significant (p < 0.001) (Figure 5).The correlations between epiphytic lichen diversity and dry branch thickness were significant (p < 0.05) (Figure 5).
A linear regression analysis was conducted between lichen species richness and single tree properties, like Dbh, tree height, crown height, and bark pH (Table 4).The models were significant but the R 2 values were low; the highest R 2 was 0.0577.There was no significant relationship between lichen richness and dry branch thickness and the total nitrogen content of bark (Table 4).
richness and individual tree properties, including the Dbh, dry branch thickness, crown height, total nitrogen, and bark pH (Figure 5).Correlations between lichen species richness and Dbh and total nitrogen of bark were not significant (Figure 5).We detected positive correlations between the epiphytic lichen richness and crown height and tree height; and negative correlations with bark pH were significant (p < 0.001) (Figure 5).The correlations between epiphytic lichen diversity and dry branch thickness were significant (p < 0.05) (Figure 5).A linear regression analysis was conducted between lichen species richness and single tree properties, like Dbh, tree height, crown height, and bark pH (Table 4).The models were significant but the R 2 values were low; the highest R 2 was 0.0577.There was no significant relationship between lichen richness and dry branch thickness and the total nitrogen content of bark (Table 4).

Epiphytic Lichen Species in Black Pine Forests
In this study, we investigated the effects of environmental factors on lichen species diversity and species composition.In the black pine forests of the Alaçam mountains, 33 epiphytic lichen species were identified.At Mt. Olympus, which shares the same latitude as the Alaçam mountains, 23 epiphytic lichen species were found [52], and 73 species were detected at Mt. Trapezitsa [53].The higher number of lichen species found at Mt. Trapezitsa compared to the other two locations could be attributed to the mixture of fir and beech species with black pines in the forests.At Mt. Uluda g, which is one of the closest rangelands to the Alaçam mountains, 20 lichen species were found [54].The low number of lichen species found in the Alaçam mountains may be due to factors related to the site conditions and the tree morphology of black pine forests.The studies conducted in mixed tree species in pine forests support this hypothesis [38,62,77].Some tree species host a greater lichen richness than black pines.For example, in a study in the Bolu mountains (Turkey), 72 lichen species were found on only five trees [78].In a study conducted in Spain, higher lichen species diversity was recorded [48].The lichen species diversity is reinforced by the exposure of tree stems to sunlight [79,80].The findings from previous studies positively support our results in the black pine forests with a mono-species composition.
An obvious summer drought occurred in the Alaçam mountains due to their long distance from the oceanic humidity-conveying winds.Black pine forests are generally found in cold and dry sites [81,82].The presence of typical crustose species in black pine forests is an indicator of these climatic conditions.In the current study, the percentage of crustose species was 27% of the total detected lichen species, which is relatively high in a particularly protected black pine forest.The branch structure affects the humidity level of the stem.Due to the sympodial structure of oak trees, the rainwater leaks through the stem, creating more favorable conditions for lichens.For pine trees, rainwater directly reaches the ground without flowing along tree branches, which in turn, produces dry substrate conditions for lichen species [77].
Hypogymnia physodes and Pseudevernia furfuracea were detected as the most widespread species at Mt. Olympus (Greece) [52].The same species showed a wide distribution in the northwestern black pine forests at Mt. Uluda g [54].Our research revealed that the most commonly distributed species in the Alaçam mountains were Pseudevernia furfuracea, Hypogymnia farinacea, Hypogymnia tubulosa, Bryoria capillaris, Platismatia glauca, Hypogymnia physodes, Ochrolechia turneri, Buellia griseovirens, Hypocenomyce scalaris, Usnea subfloridana, Evernia prunastri, Parmelia saxatilis, and Tuckermanopsis chlorophylla.P. furfuracea is a common species in black pine forests, which we also detected in the Alaçam black pine forests.The genus Hypogymnia was represented by three species in the Alaçam black pine forests, which has not been reported in previous studies.
Hypocenomyce scalaris is a species living in acidic tree bark [58].H. scalaris has been found in the oldest stand life stage (>120 years).H. scalaris has also been found on higher than 75 cm height cedar trees in Mediterranean forests [83].Sevgi et al. conducted a study in the same region and found that black pine bark pH is inversely correlated with age [84].Marmor and Randlane concluded that H. scalaris occurs more frequently on tree bark with an air-pollution-induced low pH [85].In the current study, the presence of H. scalaris in old stands revealed that H. scalaris is an indicator of this tree stand life stage.
Calicioid lichens are an indicator species of the age classes of trees.They are more common in old trees than young trees and are also known as indicators of ecological continuity [34,36].The calicioid species were assessed to indicate species diversity in pine forests [86].In the current study, Calicium glaucellum was found in 73-, 75-, and 82-year-old stands, and Chaenotheca chrysocephala was found in 78-and 82-year-old mature age class stands.The average Dbhs of the stands where C. glaucellum was found were 20, 22, and 31 cm, respectively; whereas the average Dbhs of stands where C. chrysocephala was found were 31 and 32 cm, respectively.The presence of these calicioid species indicates that age class diversity was preserved in the stand, demonstrating that old trees are protected by forest management in the black pine forests and this protection has a positive effect on the composition of epiphytic lichen species.

Stand-Level Lichen Species Composition
In studies conducted in the same region, annual precipitation was shown to affect the lichen microbiota [29,33,87].Bioclimatic variables have been used in models to explain the lichen species distribution [88].In the Mediterranean region, precipitation and temperature are prominent variables that affect the lichen distribution [89].In our study, bioclimatic variables were responsible for 0.32 of the lichen species distribution variation (Figure 2D).The results of our research revealed that Bio15, AridxThorn, and MonthTemp10 affected the epiphytic lichen species composition.Due to the summer drought in the Alaçam mountains, the humidity regime of the forest was affected by precipitation seasonality (coefficient of variation) and the index of the degree of water deficit below the water need.
We inferred that the factors affecting the humidity regime of the forest also affected the epiphytic lichen species composition.
Epiphytic lichen species diversity is affected by microclimatic conditions (air humidity, temperature, and light), and structural factors (canopy closure, vertical structure of the canopy, and shrub layer), both of which affect the microclimate in the stand [80].The topographic variables, such as Elev, CtchSlp, DwnCurv, HI, TnCurv, and Txtr, were the most important factors that explained the composition of epiphytic lichen species on the stand scale.These variables directly affected the lichen species composition due to their influence on the humidity of the site.Two important macroclimatic factors, temperature and humidity, were clearly identified among the factors influencing the composition of the epiphytic lichen community [90].Increased light penetration increases the temperature in the stand, altering the epiphytic lichen composition due to the lower humidity [91].However, Dymytrova et al. reported that the most important factors affecting the diversity of epiphytic lichens in beech forests are aspect and slope [92].Ardelean et al. noted the importance of northness [45].In our study, slope, eastness, and northness did not have any significant effect on the lichen species composition.
In our study, altitude (R = 0.26) had the strongest correlation with lichen species composition.Altitude has been the focus of many studies on the effects on lichen species composition.Dymytrova et al. stated that in beech forests (Ukraine), altitude has one of the greatest effects on lichen species density, and at higher altitudes, lichen species density increased [92].Aragon et al. [17] and Nascimbene et al. [93] reported the same result in oak forests in Spain and spruce alpine forests in Italy, respectively.Our results are in accordance with previous studies that described the effect of altitude on lichen species composition.
The stand variables that affect the climate of the site explain the diversity of epiphytic lichen species.In our study, mDbh and mCH explained more of the lichen species composition, whereas maxH, hecTBA, and sdCH did not.Dymytrova et al. stated that the mean Dbh is a major factor in the determination of lichen species diversity and composition [48,92].Those stand variables are positively correlated with the size of the crown.The crown structure varies by tree species and significantly affects species composition by supplying various conditions for epiphytic lichens [36,77].This effect is related to light availability and direct rain protection [36].As the tree crown radius increases, the shadow effect of trees increases, which increases the epiphytic lichen species composition.Our findings are in accordance with studies that indicated that the width of the crown affects the climate in the stand [36,94].We found that an increased crown radius (mCH) supports the epiphytic lichen species composition.Our results revealed that the crown size provided shelter from light and direct rain.Thus, the tree crown size positively affected the composition of the lichen species.The shadow of the high trees in the stand diminishes the light penetration and increases the humidity in the stands [95].The maximum height of the stand had an explanation capacity for the lichen species composition of 0.18.
The species response curve analyses demonstrated that the presence of some lichen species was directly affected by altitude (Figure 4A).The analysis showed that Buellia griseovirens and Ochrolechia turneri may be present above 1100 m a.s.l. in black pine forests.Giordani and Incerti stated that B. griseovirens has a crucial relationship with a cold humid climate in a study focusing on the effects of climate on lichen species at beech forests above 1000 m a.s.l.(Italy) [89].Platismatia glauca is detected above 1100 m and its presence increases in parallel to precipitation seasonality.In spruce forests, it has been detected above 1400 m a.s.l.[93].In the current study area, H. physodes was found below 1100 m.The same species has been found above 1400 m in beech and spruce forests [89,93].The Mediterranean-type climate may be responsible for that difference.The location seems to have different effects on the presence of the same lichen species at varying altitudes.
Although Ochrolechia turneri has been detected in cedar trees with more than 75 cm height in Mediterranean forests [83], in the current study area, the possibility of O. turneri presence decreased in trees with a mDbh greater than 30 cm (Figure 4B).The bark properties of varying trees may differ at a given diameter.The diameter of any tree species may vary regarding to the tree species.

Tree-Level Lichen Richness
Properties of trees, such as bark pH, structure, and branch density affect the diversity of epiphytic lichens [29,[37][38][39].The lichen species richness was lower on black pine trees than on beech trees, whereas it was higher at breast height than at the base of the stem in beech trees in pine-beech mixed forest [50].In the Mediterranean region, the number of lichen species changes due the diameter classes of cedar trees [83].In our study, there was no significant relationship between lichen species diversity and Dbh (Figure 5), whereas in some studies, significant relationships were detected between lichen species diversity and diameter breast height [34,92,[96][97][98][99][100] in forests other than black pine forests.Thus, the properties of black pine trees provide different conditions as a substrate for epiphytic lichen richness.
Due to their bark characteristics, coniferous forests generally form an acidic substrate environment for epiphytic lichen species [91] As tree diameter increases in black pine trees, the acid value decreases [79].Hauck stated that the bark pH of various conifer tree species varies between 3.0 and 4.5 [91].The bark pH of black pine in the current study area ranged between 3.2 and 4.8 [84].The average bark pH of black pine in our study was 4.0, varying between 3.3 and 4.9 (Figure 5).The linear regression between the bark pH and lichen richness had an R 2 value of 0.05 (p < 0.0001) (Table 4).The linear regression between the bark pH and the lichen richness had an R 2 value of 0.05 (p < 0.0001) (Table 4).In the study by Sevgi et. al. (2016), they found a negative correlation between bark pH and tree age.They attributed this relation to the loss of bark leaves over a certain age of black pine trees [84].Therefore, the decrease in the bark pH as the Dbh increased had a negative effect on the lichen species diversity.
The age and size of the tree affect the epiphytic lichen species diversity.There is greater epiphytic lichen species diversity on larger and older trees than in young trees [34,35,38,92].There is a significant correlation between tree height and epiphytic lichen species diversity [36].The results of our study support this relationship.
According to Ódor et al. [38], tree size does not affect epiphytic lichen species diversity.In that, this situation was reported due to the lack of old and large trees in the study area because of forest maintenance applications [38].However, in our study, the presence of trees of different ages, as a result of forest management practices in stands, enriched the composition of epiphytic lichen species.

Conclusions
Both presence-absence and frequency data were used for lichen species composition analyses.We noticed that the frequency of presence data was more powerful for the prediction of lichen species composition than presence-absence data.Topographical variations, such as elevation, catchment slope, downslope curvature, hill index, tangential curvature, and texture, were more explanatory of lichen species composition compared to bioclimatic, stand, and parent rock variables at the stand level.The bioclimatic variables Bio15, MonthTemp10, and AridxThorn and the stand variables mDbh, mCH, hecTBA, maxH, and sdCH were shown to have significant effects on lichen species composition.The lichen Bryoria implexa was found to be an indicator species for the Simav region, and Hypocenomyce scalaris was found to be an indicator species for stands older than 120 years.There was a positive relationship between tree height and lichen species composition, a negative relationship with the pH of tree bark, and no relationship with age.We suggest that the effects of topographical indices on lichen species composition should be determined within various site conditions since they have the potential to explain lichen species composition.

Figure 1 .
Figure 1.Location of the three study regions (1: Dursunbey, 2: Sındırgı, 3: Simav) of Pinus nigra forest stands in the Alaçam Mountains (top).The richness of each sampling plot is represented by the sizes and colors of circles, which represent different study regions (bottom).

Figure 1 .
Figure 1.Location of the three study regions (1: Dursunbey, 2: Sındırgı, 3: Simav) of Pinus nigra forest stands in the Alaçam Mountains (top).The richness of each sampling plot is represented by the sizes and colors of circles, which represent different study regions (bottom).

Figure 3 .
Figure 3. Ordination biplot of the stand level redundancy analysis of epiphytic lichen species.Distribution related to environmental variables: (A) LPA, (B) LFQ, (C) GLPA, and (D) GLFQ.Species are indicated by red labels and abbreviated with cepnames.Explanatory variables are indicated using blue arrows.

Figure 3 .
Figure 3. Ordination biplot of the stand level redundancy analysis of epiphytic lichen species.Distribution related to environmental variables: (A) LPA, (B) LFQ, (C) GLPA, and (D) GLFQ.Species are indicated by red labels and abbreviated with cepnames.Explanatory variables are indicated using blue arrows.Diversity 2019, 11, x FOR PEER REVIEW 11 of 23

Figure 5 .
Figure 5. Pearson pairwise correlation analysis between epiphytic lichen species diversity and tree properties.

Figure 5 .
Figure 5. Pearson pairwise correlation analysis between epiphytic lichen species diversity and tree properties.

Table 1 .
Environmental variable groups used in the study area.

Table 2 .
Results of the indicator epiphytic lichen species analysis for region and age classes.
#A refers to the probability of the suitability of studied site to the target site group, where the concerned species is present.##B refers to probability of presence of the concerned species in the target site group.Significance codes: *** p < 0.001, ** p < 0.01, and * p < 0.05.LPA: presence-absence data; LFQ: frequency data.

Table 3 .
The variance of the significant explanatory variables used in the stand-level redundancy analysis (RDA) with LPA, GLPA, LFQ, and GLFQ data sets.