Modeling of population dynamics of the protected lichen Lobaria pulmonaria (L.) Hoffm. in the Caucasus

Reduction of biological diversity of lichens, reduction of the distribution of rare species and their disappearance due to habitat disturbance are significant problems in the Caucasus. The aim was to study the main patterns of distribution of the rare lichen species Lobaria pulmonaria (L.) Hoffm. and included identification the main abiotic factors affecting the distribution of the species in region. We modeled the current habitats of Lobaria pulmonaria in the Caucasus by using the Maxent method (Maximum entropy modelling). The most suitable for distribution of the lichen were the wettest areas of the southern macroslope of the Greater Caucasus. The center of the predicted range of L. pulmonaria was currently located on the Black Sea coast, in the mid-mountain areas of Krasnodar region, Georgia and Abkhazia. The minimum probability of finding the species predicted in relatively arid areas with a more continental climate in the Central and, especially, Eastern Caucasus. Temperature and orographic (Topographical Ruggedness Index, topographical humidity Index) factors are also important in the distribution of the studied species.


Introduction
The lichen Lobaria pulmonaria (L.) Hoffm. is a species with a wide range, but currently protected in Russia [1], as well as in a number of European countries and Canada [2][3][4]. The species grows in boreal and temperate zones, non-tropical rainforests, as well as in mountainous areas and areas with a maritime climate. Many scientists noted that the lichen range is currently decreasing a result of deforestation and air pollution, but we do not know the degree of the species vulnerability to these influences. In the Caucasus, L. pulmonaria occurs mainly on trunks of the trees in forest regions, but the cenotic confinement of L. pulmonaria poorly studied. It known that the lichen grows in different types of forest communities: spruce forests, old growth aspen forests, coniferous and smallleaved forests, as well as in forests with broad-leaved species [5]. The study of lichen cenotic confinement is important for monitoring the status of this rare species, searching for new habitats, and developing strategies for its protection. But it's also important to study the spatial distribution of L. pulmonaria and to predict its possible locations the Caucasus.
The aim was to study the main patterns of L. pulmonaria distribution in the Caucasus including the Caucasian floristic province of the Southern European mountain forests region. The objectives were a) to identify the main abiotic factors that influence the distribution of the lichen in the study area by using the Maximum entropy modelling methods; b) to receive a map of the current potential distribution of L. pulmonaria in the Caucasus.

Study area
The field of research covered the Caucasian mountain country in the following geographical boundaries: in the east, along the coast of the Caspian Sea; in the north, from the Caspian Sea along the Kumo-Manych Depression and further along the northern border of Krasnodar krai to the Black Sea; in the west, along the Black Sea coast; and in the south, by the borders of Georgia, Armenia, and Azerbaijan with Turkey and Iran. The total area of the analyzed territory amounted to 390 thous. km 2 ( Figure 1).
The analysis of the L. pulmonaria was carried out within the Kubanskiy (Western Caucasus) and Elbrusskiy (most of the Central Caucasus), Terskiy (Eastern Caucasus) variants of vertical zonation of the Western North Caucasus zonality type and Eastern North Caucasus zonality type, respectively [6]. The territory of Western North Caucasus zonality type stretches from the eastern shore of the Azov Sea to the Teberdin-Daut upland and the Stavropol Upland Ridge, which serves as a watershed between the Black Sea -Azov and Caspian basins. The gradient belt of this zonality type is a barriersectoral variation of the steppe belt type, influenced by the littoral foothills of the Caspian basin [6]. The altitudinal zonality spectrum developed here under the influence of the steppe latitudinal zone. There is also an increased wetting pattern associated with the influence of Black Sea -Mediterranean cyclones and the proximity of the Atlantic Ocean, determines many of the ecological features that determine the existence of distinctive flora and fauna found here. The Stavropol Upland, being a barrier to westerly moist air currents, contributes to the considerable precipitation in this area. This creates conditions characteristic of for pre-oceanic sectors, and therefore the type features a weakly continental "littoral" subtype that has one, the Kubanskiy variant of vertical zonation. In the foothills, the Kubanskiy variant corresponds to West Kuban sloping alluvial plain. The mountainous part represented by the western spur of the Caucasus Mountain range comprising high-mountain and medium-mountain ridges.
Western Caucasus characterized by humid climate due to its location in a zone of Mediterranean and Atlantic cyclones. The Stavropol Upland, which acts as a barrier to them, contributes to the "barrier foothill" effect. The amount of precipitation on the Azov-Kuban lowland varies in precipitation between 500-600 mm and 700-800 mm per year. In the mountains, it increases sharply. In particular, in the Arkhyz area the amount of precipitation reaches 863 mm. Winter temperatures in the middle reaches of the Kuban river in January are, on average, -2…+3°C. Spring begins in the second or third decade of February. Summers are also much milder here than in the eastern half of the Caucasus. The average daily temperature in July is +22.5…+23.5°C. Towards the sea, there is a gradual increase in precipitation.
Characteristic features of Kubanskiy variant of vertical zonation are 1) the presence of thick, high-tree mesophilic forest ecosystems occupying foothills and mid-mountains and have an intensive reforestation process; 2) widespread development of normal subalpine meadows and subalpine mixed grasses (relict systems of the Tertiary period); 3) almost complete absence of mountain steppes.
Elbrusskiy variant of vertical zonation is located in the Pyatigorie-Elbrus area. Its north-western boundary coincides with the type boundary, while its southeastern boundary runs along Dykh-Tau -Karakaya -lower reaches of Baksan River. The main feature of the variant of vertical zonation is that the Greater Caucasus mountain range blocks the influence of moist Mediterranean -Black Sea air currents, while the advanced ridges of the Rocky and Cretaceous gorges block the influence of moist air currents. Cretaceous mountain ranges are gently sloping and do not prevent the Caspian lowland's dry winds from freely penetrating the mountains. Hence the sharp xerophitization of all the southern landscapes and the dryness and continentally of the local climate in the mountains. The peculiarities of composition and distribution of plant species are a response to the permanent impact of semidesert conditions of the Caspian Lowland and the highly developed glaciation of the highlands of the Central Caucasus. Hence the wide steppification of the subalpine belt and the alpine belt in some places, as well as the loss of the belt of dark coniferous forests and the belt of broad-leaved forests.
Terskiy variant of the vertical zonation is located in the Terek and Argun River basins. In the northwest the boundary of the variant lies along Dykh-Tau -Karakaya -lower reaches of Baksan River, and in the south -up to the crest of Andi Range, and the Greater Caucasus Range. The Terek variant characterized by a strong convergence of the Main, Lateral, Rocky and Cretaceous Ranges, all with a sharp profile, high and sharp peaks, with steep, stony and rugged slopes. Increasing heights of the advanced ridges as well as Terskiy and Sunzhenskiy Ridges sharply increase their role as a barrier on route of dry winds from Caspian Sea that considerably weakens influence of north-eastern dry winds on mountainous landscapes. This, in turn, enables mountain physical-geographical altitudinal patterns to manifest themselves more freely. The peculiarity of the natural and climatic conditions of this region determines the wide distribution here, as in the Kubanskiy variant, of broad-leaved forests, however, more arid and depleted of tertiary elements. Moreover, the belt of dark coniferous forests typical of the Kubanskiy variant is absent. Occasionally on the more arid slopes isolated islets of pine forests. In connection with the above, the species composition of the Caucasian mesophilic groups is again expanding, but unlike the Kubanskiy variant, they are not numerous and their distribution and choice of biotopes is limited. There is an increase in the influence of xerophilic, steppe and semi-desert species, an increase in their upper habitat limit, although not to the same extent as in Elbrusskiy variant.

Data collection and measurements
During the 2008-2020 expeditions, we identified 48 presence points of L. pulmonaria in the Caucasus. Each habitat described using several indicators: GPS coordinates (accuracy 3 m), height above sea level, ecological substrate. The location of each sample on the substrate (height and exposure) recorded. The geographical coordinates of the localities were determined using a GPS navigator "Garmin".

Assessment and analysis of environmental predictors
We used an extended set of ENVIREM (ENVIronmental Rasters for Ecological Modeling) [7] bioclimatic and topographic modeling layers to estimate the main factors driving the distribution of L. pulmonaria, largely complementing the widely used WorldClim database. The set consists of 2 topographical variables and 16 bioclimatic variables, including information about evaporation processes from the soil surface and vegetation [8]. The connections between the ENVIREM layers and the physiological processes in plants determines the feasibility of these predictors use for modeling the areas of plant species.
At the first stage of the assessment and analysis of environmental predictors, we used full ENVIREM set of bioclimatic and topographic variables.
At the second stage, the VIF test (Variance Inflation Factor) in the R program was used to assess the collinearity between the initial environmental predictors when modeling species distribution, which, among other things, the presence of latent correlation. The values obtained for a VIF threshold ≤ 10 allowed us to identify eight predictors for subsequent modeling of the spatial distribution of L. pulmonaria (Table 1).
At the third stage, for comparative analysis, we conducted a factor analysis (PCA) of the ENVIREM environmental predictors. PCA also allows you to evaluate and eliminate the collinearity between the ENVIREM variables. PCA revealed three main components with eigenvalues of factors greater than 1 ( Table 2). These eigenvalues of the factors (the first three orthogonal axes of PCA) later used as layers for modeling.

Maxent modeling and model evaluation
Maxent (Maxent software for species habitat modeling) [9,10] was used to model the spatial distribution of L. pulmonaria within the study area.  To obtain an adequate model, the calculation was performed by five repetitions (to obtain a test sample of 20-25%), automatic functions (Auto LQP), using 500 iterations for each pixel of the entire analyzed area. The quality of the "candidate models" was evaluated by the values of the Akaike information criterion for small samples AICc (Akaike's information criterion corrected) and the Bayesian information criterion (Bayesian information criterion), aimed at selecting models based on a compromise between their complexity and accuracy. In addition, we used AUC (Area under the curve), which indicates the results of statistical analysis of the coincidence of models constructed from test and training data (evaluation of the model accuracy). The final (optimal) models in each case selected according to the minimum values of the AICc and BIC at high AUC values (Table 3).  Based on AUC, AICc, and BIC values, all models had adequate values. The AUCs showed good values for the training and test models. However, of the three data sets presented, the ALL_ENVIREM model data set had the highest AUC, while the full set also had the lowest information criteria values (Figure 2).

Fig. 2. AUC, AICc, and BIC values depending on the model type
If we proceed from the accepted criteria for model selection, then in our case the most adequate choice between model complexity and accuracy was in the area of minimum values of information criteria and maximum values of AUC. In our models, ALL_ENVIREM, model 4 met such criteria. However, in such a case, accounting for correlated variables might cause overfitting of the model. Therefore, in this situation, it was better objectively use the axes of factor analysis, which eliminated the problem of collinearity and allowed to calculate models based on the weights of factors in the cells of the raster. Proceeding from the above, the most adequate was PCA_ENVIREM, model 2, for which the description of the spatial distribution given.

Results
As shown above, the main variability of environmental parameters (92.7% of cumulative variation) at the presence points of L. pulmonaria in the Caucasus explained by the first three PCA axes. The analysis of factor loadings showed that the first main axis (about 62% of variables variance) is mainly formed by seven predictors (Table 2). Two highly correlated variables of the first axis characterized potential evapotranspiration (an indicator of the maximum amount of moisture evaporated by plants from the reference surface in the absence of moisture deficiency): annualPET is the average annual potentional evapotranspiration, and PETColdestQuarter is the average monthly potentional evapotranspiration for the coldest quarter. Evapotranspiration depends on the solar radiation, air temperature and wind speed, which promote plant transpiration, as well as on the frequency of precipitation in the region [11]. Five variables are directly connected to the temperature of the surface layer of the atmosphere: growingDegDays0 и growingDegDays5 are the sum of the average monthly temperature for the months with an average temperature above 0 and 5°C, respectively, multiplied by the number of days; maxTempWarmestMonth is the maximum temperature of the warmest month; monthCountByTemp10 -the number of months during which the average air temperature exceeds 10°C; thermInd is the index of compensated thermality (depends on the sum of the average annual temperature range, minimum and maximum temperatures of the coldest month).
The main variables of the second PCA axis characterized the water regime of the study area: aridityIndexThornthwaite is the Thornthwaite aridity index of climate; climaticMoistureIndex is the index of relative humidity and aridity of the climate; embergerQ is the Emberger pluviothermic coefficient, developed for the differentiation of the Mediterranean climate type (associated with the annual potential evapotranspiration and extreme annual thermal amplitude).
The third PCA axis combined orographic parameters: TRI (Topographic Ruggedness Index) -an index of terrain roughness, which characterizes the local vertical dissection of the terrain [12]; topoWet -a topographic humidity index used to study the moisture level of substrates in different parts of the slopes in the mountains [13].
The map constructed based on the final Maxent model (PCA_ENVIREM, model 2) showed that under modern climatic conditions, L. pulmonaria had a wide potential distribution in the Caucasus (Figure 3). The areas that are optimal for the growth of the species were mainly concentrated on the southern macroslope of the Great Caucasus. The probability of L. pulmonaria findings was higher in relatively wetter mid-mountain areas, noticeably decreasing to the east of the central part of the Greater Caucasus. The most extensive areas of optimal habitats of the species covered mountainous regions of Krasnodar region, Georgia and Abkhazia. On the Black Sea coast of the Caucasus, the area of optimal habitats for the species decreased in the direction from the southeast (areas with a humid subtropical climate) to the northwest (areas with a more temperate climate). On the northern macroslope, L. pulmonaria had a high potential for distribution from the mid-mountains to the highlands in the western part of the Greater Caucasus. Areas with a species finding probability above 80% were concentrated in Krasnodar region, Republic of Adygea, Karachai-Cherkess Republic. Rarely, the areas optimal for the species growth found in the relatively arid continental conditions of the Central and, especially, the Eastern Caucasus, including the coast of the Caspian Sea.
In the central and eastern parts of both macroslopes of the Greater Caucasus, optimal habitats for L. pulmonaria were limited to small areas of the wettest mountain areas, the area of which decreases from west to east along the gradient of increasing climate aridity. The high mountainous areas of the predicted range of L. pulmonaria in the Caucasus were concentrated around the valleys of large rivers. 0.5-0.8 and above 0.8, the probability for suitable and optimal habitats, respectively The significance of the water regime in the distribution of L. pulmonaria in the Caucasus confirmed by the analysis of the contribution of PCA factors to the construction of the prognostic Maxent model of the species distribution ( Table 5). The largest percentage contribution to the construction of the final Maxent model (PCA_ENVIREM, model 2), both with an independent influence (percentage contribution) and considering the correlation with other factors (permutation impotance), was made by Factor 2. It combined, as shown above, the parameters of the water regime -the Thornthwaite aridity index, the index of relative humidity and aridity of climate, and the Emberger pluviothermal coefficient. At the same time, the ariditylndexThornthwaite and continentality characterized by a negative value of factor loads, while the embergerQ and climaticmoistureldex had a positive correlation with Factor 2. Accordingly, Factor 2 characterized as a "factor of the moisture availability" of the environment in the predicted habitats of L. pulmonaria in the Caucasus. This largely explained the potential distribution of the species in the Caucasus with a pronounced confinement to the wettest territories, the area of which reduced in the direction from west to east on the gradient of increasing climate aridity.
The second most important complex Factor 3, used as a predictor in modeling, included both topographic parameters. So, Factor 3 characterized as a "topographic factor". The terrain roughness index TRI, which increases with increasing slope steepness, characterized by a negative factor loading, while the topographic humidity index topoWet positively correlated with Factor 3. Accordingly, the most suitable places for L. pulmonaria in mountainous areas were non-steep slopes with moist soils.
"Temperature" Factor 1, despite explaining a most of environmental variability at the presence points of L. pulmonaria, was of the least importance in predicting the spatial distribution of the species in the Caucasus.

Discussion
The purpose of this study was to model the spatial distribution of L. pulmonaria and to identify the main abiotic factors affecting the distribution of the species in the Central Caucasus. We revealed the significant effect of precipitation and temperature on the distribution of L. pulmonaria. These results are in line with previous reports about the ecology of the studied species, according to which L. pulmonaria tends to wet forest areas. For example, Eaton and Christopher [14] used the generalised linear mixed models to compare thallus growth to a suite of climatic variables (WorldClim dataset). They found a significant relationship between thallus growth of L. pulmonaria in forest habitats and macroclimatic variables, among which the main ones were total precipitation and annual mean temperature [14]. Muir et al. [15] concluded that growth in biomass of L. pulmonaria depended on humidity and varied seasonally in areas with large seasonal variations in precipitation. Pannewitz et al. [16] also reported that alterations in the temperature and water regimes are especially detrimental to this shade adapted lichen.
Modeling the potential distribution of lichens of the Xanthoparmelia pulla group (Xanthoparmelia delisei, X. loxodes, and X. verruculifera) in Central Europe also showed that potential distribution of the species strongly influenced by precipitation-related variables [17]. Our results are in line with K. Pearson et al. [18], which observed that, according to Maximum entropy modelling, the predicted distribution of lichen Fuscopannaria leucosticta in Canada most affected by temperature. On the other hand, E. Holt et al. [19] concluded that lichens influenced by local factors (e.g., disturbance, competition) which may mask any largescale patterning such as precipitation, topography. In general, as Ch.J. Ellis [20] showed, there exists a strong framework for using bioclimatic models to better understand lichen distributions, but lichenologists should seek to ensure that model of species distribution must be supported by prior knowledge of lichen functional ecology. The author also concluded that larger-scale climate may interact with small-scale variables relating to topography, forest stand structure, etc. [15].
According to D. Stoykov [21] and Jüriado and Jaan [22], L. pulmonaria is an indicator of old growth forests, mainly found on the trunks or bark of mature broadleaved trees and less often on conifers (Abies sp., Picea sp. and Pinus sp.). It prefers shady habitats in undisturbed forests [16]. On plains, the species prefers ridges and often grows in floodplain forests, on alder and beech, and in mountains, its habitat is beech forests [23,24]. Since lichens must "follow their phorophytes" to survive [25], the obtained map of the potential distribution of L. pulmonaria in the Caucasus largely corresponds to the map of the distribution of wet oldgrowth forests in the Caucasus. Thus, the preservation of undisturbed old-growth forests in the Caucasus is very important for the conservation of rare lichen species such as L. pulmonaria.

Conclusions
The analysis of the main ecological predictors, which determined the potential distribution of L. pulmonaria in the Caucasus, confirmed the dependence of the studied species on the parameters of water and temperature regimes. At the same time, the predicted spatial distribution of the species mostly depended on the complex factor of moisture availability of the study area (included the Thornthwaite climate aridity index, the index of relative humidity and aridity of the climate, and the Emberger pluviothermal coefficient). Lobaria pulmonaria is widely represented on maps of potential distribution within the boundaries of the northern and southern macroscopes of the Greater Caucasus. The center of the current potential range of L. pulmonaria is located in the wettest regions of the Caucasus -on the Black Sea coast, and in the mid-mountain regions of Krasnodar region, Georgia and Abkhazia. The minimum probability of the species findings predicted in relatively arid areas with a more continental climate in the Central and, especially, Eastern Caucasus. Temperature and orographic parameters were also significant in the distribution of the studied species. The high-altitude areas of the predicted range of L. pulmonaria in the Caucasus are concentrated around the valleys of large rivers.
This study identified L. pulmonaria habitats and species-habitat relationships, providing conservationists with baseline data that can help locate additional populations and provide a better understanding of the ecological features of this lichen, which will help develop robust conservation strategies for this rare species.
The studies were carried out as part of state assignment no. 075-00347-19-00 on the topic "Patterns of the Spatiotemporal Dynamics of Meadow and Forest Ecosystems in Mountainous Areas (Russian Western and Central Caucasus). "