Population and Conservation Status of Buxbaumia viridis (DC.) Moug. & Nestl. in Romania

The field survey made in the last 20 years revealed that large areas of Romania’s territory are still unexplored concerning moss distribution. The long-term research in natural and mature spruce forest habitats of this country shows that many sites are well protected, a status that is also confirmed by reports of Buxbaumia viridis. However, many other sites where this species was identified still lack legal protection. We also generated a potential distribution of the species based on an ensemble model, useful in guiding extensive field surveys and also management and conservation measures. In a country where the volume of wood cut by illegal logging is larger than the legal one, according to official data, it is very important that all habitats for B. viridis be included in protected areas. Our paper brings important data to aid in this goal.


Introduction
Buxbaumia viridis (DC.) Moug. & Nestl. is a European protected species, listed in Annex II of the Habitats Directive (DIRECTIVE 92/43/CEE) [1]. This species is characteristic of natural, mature spruce forests. It grows especially on dead wood, but occasionally it can also be found on soil, provided wood remains exist [2]. It prefers wood in an advanced stage of decay, being found on large fallen spruce trunks, especially near streams, due to its need for high humidity. As a result of its ecology, B. viridis is considered an indicator species for spruce forests in good and very good conservation conditions [3] but can also be found in managed forests [2,4].
B. viridis is a Least Concern (LC) moss of Europe but is a threatened species in many countries: it is Critically Endangered (CR) in Serbia, Montenegro, and Italy, Endangered (EN) in Finland, Greece, Hungary, Romania and Vulnerable (VU) in Andorra, Great Britain, Spain, Czech Republic, Slovakia and Estonia [5].
According to the Ministry of Environment, Waters and Forests, it is estimated that about 38.6 million m 2 of wood are logged every year in Romania, but only 18.5 million are legally extracted [6]. Even if Romania has a large area of natural forest, illegal cutting represents the biggest threat to mosses such as B. viridis.
B. viridis has been known in Romania since the 19th century, and there are many confirmations of the species' presence on Romanian territory. However, recent research has shown that there are still many areas that have not yet been fully explored, and B. viridis can be found there, especially if suitable habitats exist but have not been previously detected.

Current Species Distribution
In this study, we present an updated distribution of the moss B. viridis in Romania, with recent field-collected data, first-time reports for 15 mountain regions and confirmations for other mountainous areas. Buxbaumia viridis was found in mountain massifs and lowland areas, and some mountain massifs such as Ceahlău, Vrancea, Penteleu, Siriu, Ciucaș, Grohotiș, Baiului, Postăvaru, Leaota, Făgăraș, Păpușa, Șureanu, Retezat and Țarcu represent new reports for this species. We assembled a database with 195 spatial records by combining 159 original and 36 published distribution information. We found sporophytes during snowy winter, demonstrating that it is not impossible to find them in this period [14]. All known distribution records for B. viridis in Romania are presented in the subsequent maps (Figures 1-4).         [40,41], but from the paper results, the records are from the northern part of Pop Ivan Peak, which is on the territory of Ukraine.
Also, some doubtful reports from Iezer-Păpușa Mountains [10] were excluded from this study.
New records of Buxbaumia viridis:  Doubtful records for Romania: B. viridis was reported from Pop Ivan Peak, Maramures , Mountains [40,41], but from the paper results, the records are from the northern part of Pop Ivan Peak, which is on the territory of Ukraine.
Also, some doubtful reports from Iezer-Păpus , a Mountains [10] were excluded from this study.
New records of Buxbaumia viridis:

Habitat Suitability
The distribution modeling displayed good metric performance with mean AUC values between 0.81 and 0.89 and mean TSS values in the range of 0.58 to 0.67 (Table 1). ROC plots also displayed good performances by all modeling algorithms, with AUC values > 0.8 (Figures S1-S3). The highest contributing variables were elevation (16%), precipitation of driest quarter (12.9%), and mean temperature of coldest quarter (11.4%), followed by precipitation of warmest quarter (8.1%), percentage of tree cover (3.6%), heat load index (3.3%), leaf area index (3%) and Euclidean distance to the nearest river with the lowest contribution in constructing the model (1.2%) (Table 2, Figure S4). The response curves show that the species' range encompasses mainly high-altitude areas, generally located between 1000 and 1500 m, with high levels of precipitation during the driest quarter of the year (above 150 mm), low mean temperatures during the coldest quarter (below −5 • C) and high precipitations during the warmest quarter (about 300 mm) ( Figure S4); increased tree cover, slopes with low heat load index, moderate leaf area index are also important in defining the environmental space for the green-shield moss ( Figure S5).
With regard to the known distribution of Buxbaumia viridis, the present model of habitat suitability showed remarkable congruence. The ensemble model highlighted a high density of suitable habitats for B. viridis in most of the mountain massifs where the species was found, such as Retezat, Bucegi, Făgăras , , Piatra Craiului and the western part of Apuseni Mountains ( Figure 10). Furthermore, wide areas of highly suitable habitats for this moss have been identified in Gutâi, T , ibles and Maramures , mountains in the north, but also Gurghiu and Semenic mountains. The central part of Eastern Romanian Carpathians has mostly isolated areas of suitable low to medium habitats, as do the Poiana Ruscă Mountains and the eastern part of the Apuseni Mountains. Out of all 167 physiographic units identified as holding some potential areas for the species' range by the threshold model, 129 have no known occurrence records.

OR PEER REVIEW 11 of 18
The distribution modeling displayed good metric performance with mean AUC values between 0.81 and 0.89 and mean TSS values in the range of 0.58 to 0.67 (Table 1). ROC plots also displayed good performances by all modeling algorithms, with AUC values > 0.8 ( Figures S1-S3). The highest contributing variables were elevation (16%), precipitation of driest quarter (12.9%), and mean temperature of coldest quarter (11.4%), followed by precipitation of warmest quarter (8.1%), percentage of tree cover (3.6%), heat load index (3.3%), leaf area index (3%) and Euclidean distance to the nearest river with the lowest contribution in constructing the model (1.2%) (Table 2, Figure S4). The response curves show that the species' range encompasses mainly high-altitude areas, generally located between 1000 and 1500 m, with high levels of precipitation during the driest quarter of the year (above 150 mm), low mean temperatures during the coldest quarter (below −5 °C) and high precipitations during the warmest quarter (about 300 mm) ( Figure S4); increased tree cover, slopes with low heat load index, moderate leaf area index are also important in defining the environmental space for the green-shield moss ( Figure S5).
With regard to the known distribution of Buxbaumia viridis, the present model of habitat suitability showed remarkable congruence. The ensemble model highlighted a high density of suitable habitats for B. viridis in most of the mountain massifs where the species was found, such as Retezat, Bucegi, Făgăraș, Piatra Craiului and the western part of Apuseni Mountains ( Figure 10). Furthermore, wide areas of highly suitable habitats for this moss have been identified in Gutâi, Țibles and Maramureș mountains in the north, but also Gurghiu and Semenic mountains. The central part of Eastern Romanian Carpathians has mostly isolated areas of suitable low to medium habitats, as do the Poiana Ruscă Mountains and the eastern part of the Apuseni Mountains. Out of all 167 physiographic units identified as holding some potential areas for the species' range by the threshold model, 129 have no known occurrence records. The field surveys carried out resulted in 218 microhabitats colonized by B. viridis at the national level. From these, 172 microhabitats have been recorded in this study, with a total of 1534 sporophytes. In Romania, B. viridis is characterized by the fact that it colonizes The field surveys carried out resulted in 218 microhabitats colonized by B. viridis at the national level. From these, 172 microhabitats have been recorded in this study, with a total of 1534 sporophytes. In Romania, B. viridis is characterized by the fact that it colonizes spruce deadwood [3], as opposed to its occurrence on soil in Hungary [2]. Our field activities identified 162 microhabitats (>94%) belonging to coniferous deadwood, mainly Picea abies (L.) H. Karst., five microhabitats were beech deadwood (Fagus sylvatica L.), two microhabitats were white alder deadwood (Alnus incana (L.) Moench) and three microhabitats were on soil.

Population and Conservation Status
Considering that the Area of Occupancy (AOO) according to our study is 268 km 2 , representing approximately 1% of the total potential distribution range of the species, we can conclude that the total population of Buxbaumia viridis in Romania is over 150,000 sporophytes distributed in more than 17,000 microhabitats.
We propose that the conservation status of Buxbaumia viridis in Romania should be changed from Endangered-EN A3c; C1 [5,12] to Vulnerable-VU A3c.

Discussion
Our research showed, as we expected that in places where the natural and mature spruce forest habitats are present, the probability that Buxbaumia viridis is present is very high. Moreover, the area covered by spruce sums up to 1.37 million hectares, as estimated in the latest National Forest Inventory [42]; hence a large area with potential habitat in Romania is still uncovered by research.
Although it is an ephemeral species, B. viridis appears from August-September and matures the following year, sporulating from May-July [43]; populations can be found years in a row in the same place. For example, in the Giumalau Massif, the species was identified on a spruce log in 2007 and 2008 [39], and in 2018 it was found on the same log, in the same area of the log, at the base, so after a period of 12 generations.
The compiled distribution map, overlaid with the borders of all the sites of Community importance (SCI, as part of the NATURA 2000 network), revealed that 20 more areas have no protection at all, and these will be part of the new NATURA 2000 site proposals.
According to the results from the ensemble modeling, elevation, precipitation of the driest quarter, and mean temperature of the coldest quarter serve as the best indicators of Buxbaumia viridis occurrence. We can observe that the favored high-altitude areas typically correlate with the altitudinal range of mixed and spruce forests; hence there might be a dependence on habitat type, also supported by other studies [44,45]. Using maximum entropy modeling,Číhal [44] also argues that the found dependency of B. viridis to habitat type is most likely mostly connected to the species' requirement for a sufficient quantity of decaying wood [46]. However, when focusing on desiccation, Kropik [47] found that this variable outperforms decaying wood in terms of predicting power for explaining the occurrence of this moss in Austria. This result indicates that climate can have a more substantial effect on the distribution. Our data suggest the same, as precipitation of the driest quarter and mean temperature of the coldest quarter accompanied by precipitation of the warmest quarter are the important consequent variables after elevation. Increased precipitation (in the driest and warmest areas) has a direct influence on avoiding the desiccation of not just sporophytes [45,46] but also of the less tolerant spores [48]. The reason for the species' dependency on lower average winter temperatures is uncertain, as also found and stated byČíhal [44]. Only 3 • C of winter warming in several experimental plots in the UK caused a cover increase in certain moss species and a decrease in cover in other species, under a general decrease of bryophytes species richness. Little is known about bryophytes' spore dormancy [49] or spore development that could explain this dependency. Other causes could be indirect.
In the Pyrenees, the effects of organisms grazing on sporophyte were observed in spring, sometimes with slugs from the Arion genus being responsible [50]. However, slugs can be active even in winter, and the activity of some species increases when the soil temperature is above 10 • C [51]; decaying wood is also known to act like temperature microhabitats refuges for soil invertebrates. For example, adults of Arion distinctus Mabille, 1898 are thought to survive in winter and even lay eggs throughout these months until spring [52]. Nevertheless, despite their freeze tolerance, slug species from the same genus do not survive below −3 • C [53]. Therefore, we could argue that low mean temperatures below −5 • C might reduce the slugs' survival and density, hence the grazing intensity on B. viridis sporophytes. Both spore dormancy/development and bryophagy effect over this moss during winter, in correlation with low temperatures, are worth further exploring.
The next three variables in order of importance in producing the model are the percentage of tree cover, heat load index, and leaf area index. Slopes with low heat load index also indicate shadier and more humid forest floor conditions, which is in agreement with findings about "Northness" in other papers [45,54], but with a finer prediction due to the calculation method [55]. Spitale and Mair [54] also found that this moss species prefers closed-canopy forests. Canopy closure, percentage tree cover and leaf area index are structural variables describing the canopy. The canopy has an essential impact on understory composition and species cover in a forest [56] by influencing light, humidity and even soil pH. Our data showed that this moss is most likely found in areas with increased tree cover but a moderate leaf area index. In our field experience, coniferous forests with higher leaf area index tend to accumulate a thicker layer of needles that seemed to hinder gametophyte development.
As we suspected, the model for habitat suitability detected physiographic units (mountains and lowland areas) within the species range, some with high suitability, such as Gurghiu and Semenic. These areas, as well as the northern part of the Romanian Carpathians (Maramures , mountains), should be the first target in future field studies. The record from Stânis , oara Mountains, Găines , ti forest chalet is the only one that is situated outside the species' range. With new field data, especially from the central outer part of the Eastern Romanian Carpathians, the model results could be improved, but it is also worth confirming the record again, as the last mention was from four decades ago [23]; Previous research indicated that deadwood amount is an important variable influencing B. viridis distribution predictions [46,47,54], which could improve distribution modeling. However, spatial databases and statistics regarding forest deadwood were only recently published, with a low resolution of 16 km, at the European level and did not cover Romanian territory [57]. Furthermore, some anthropogenic-related disturbances like air pollution were found to have negligible influence on this moss distribution [44]. Other categorial predictors like land use (Corine Land Cover) or forest management methods are sometimes found significant [44,45], but these also determine deadwood volume, at least in some forest types [58]. In addition, land use data are usually highly correlated with other environmental predictors, such as the climatic envelope. Moreover, categorical predictors are difficult to include in species distribution modeling since not all algorithms successfully manage them.
The new reports of B. viridis from areas not currently protected, such as Capăt , ânii, Baiului, Grohotişului and Ous , or Mountains, offer support for the designation of new Natura 2000 sites, as well as the opportunity to expand the limits of existing Natura 2000 sites of community importance, especially ROSCI0101 Larion, ROSCI0019 Călimani-Gurghiu, ROSCI0208 Putna-Vrancea, ROSCI0190 Penteleu, ROSCI0229 Siriu and ROSCI0207 Postăvarul.

Occurrence Records
An extensive field survey for B. viridis was conducted by the authors in South-Eastern Carpathians in the last 20 years; the following data were collected for each occurrence: geographic coordinates (decimal degrees) and elevation using a GPS receiver, the number of microhabitats (logs or 1 m 2 area for soil) and number of sporophytes in all identified presence locations. A critical review of scientific publications was carried out, and distribution data were registered in a database. When no coordinates were provided, toponyms were georeferenced using Google Earth Pro version 7.3.6.9345, Google Maps version 2022, and military survey maps for Romania with an assumed error of less than 1 km. Distribution maps were created using the collected spatial data in ArcGIS 10.7.1 [59].

Model Preparation and Procedure
We used both original and published distribution records for the green-shield moss (B. viridis); before modeling, we rarefied the occurrence dataset using SDMToolbox v2.5 add-on [60] in ArcGIS 10.7.1 and a distance of 1 km to match the spatial resolution of the input predictors; the resulting dataset used for modeling comprised 95 occurrence records.
Regarding environmental predictors, we started with a package that contained 12 variables, describing climatic, ecological and geomorphological conditions in the area where the species is distributed (Table 2). For the bioclimatic predictor, we used the secondgeneration baseline (v2.1) [61], representative of the 1970-2000 period. Heat Load Index (HLI) was calculated using the Geomorphometry and Gradient Metric Toolbox v2.0-0 [62]. All variables were used at a resolution of 30" (~1 km).
We first split the environmental rasters into two sets of predictors: (1) one set for training, clipped to the minimum convex polygon enclosing the input features (convex hull) plus a generic buffer of 10 km, and (2) a second set for predictors, clipped to the border of Romania plus a generic buffer of 10 km. This phase was carried out in ArcGIS 10.7.1 using the SDMToolbox v2.5 add-on.
Before the modeling phase, multicollinearity among variables was tested using the Variance Inflation Factor (VIF) implemented in the usdm package [63] in R version 4.2.1 [64] and excluded variables with VIF > 10 [65]. The final list of predictors used is presented in Table 2.
We employed ensemble modeling in R using the sdm package [66] to predict habitat suitability for green-shield moss (B. viridis) in Romania. Models were generated using the training predictors and six methods, GLM (Generalized Linear Models), RF (Random Forests), Maxent, BRT (Boosted Regression Trees), MARS (Multivariate Adaptive Regression Spline) and SVM (Support Vector Machine), and then projected onto the entire environmental space (as defined above). Settings included 10.000 random background points, 5-fold cross-validation with 10 repetitions, and a test percentage of 30%.
In total, 300 models were generated and evaluated based on AUC and TSS. The ensemble was generated from models with AUC values greater than the mean plus half a standard deviation or models with TSS values greater than the mean plus half a standard deviation [67]. A threshold model was also created based on the TSS cut-off value.

Conservation Status Assessment
Conservation status was reevaluated based on the new distribution records from the current article and the IUCN methodology. AOO was calculated in ArcGIS 10.7.1 using a 2 km square grid, as recommended in the IUCN Red List Categories and Criteria [68].

Conclusions
This study can be used to designate new NATURA 2000 sites or enlarge the areas of the actual sites for the B. viridis species in Romania, which would contribute to the fulfillment of the target assumed by The European Commission through the EU Biodiversity Strategy for 2030: Bringing nature back into our lives in 2020, to establish a protected area for at least 30% of land in Europe [69].
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12030473/s1, Figure S1: ROC-AUC curves for the GLM and RF models; Figure S2: ROC-AUC curves for the Maxent and BRT models; Figure S3: ROC-AUC curves for the MARS and SVM models; Figure S4: Relative variable importance for all models across all modeling methods (SD bars included); Figure S5