Application of Phytosociological Information in the Evaluation of the Management of Protected Areas

The classification system of plant communities using phytosociological methods can be applied to their conservation in protected areas, as well as in establishing adequate protections and granting legal status to such areas. A new integrative index is developed to classify plant communities for the evaluation of the conservation status of protected areas, obtained from the product of three statistical indices of diversity: Syntaxonomic Distinctness, Rarefaction and Areas Prioritisation, which has been named DRA (acronym of the three indices used). The DRA is used to assess whether the status granted to Protected Areas matches the values provided by the plant communities within them and which were the basis for the identification and description of the Habitats of Community Interest (Habitats Directive—92/43/CEE). The proposed method was applied to the network of protected natural areas on the Andalusian coast, including 14 areas with different protection status, where, once the plant communities they contain were identified, the DRA index was applied to each of them and compared with the Legal Protection Index, i.e., the current protection regime; it becomes clear, objectively, that not all the statuses assigned, whether the IUCN criteria or those of the Andalusian government, correspond to the real levels of protection they should have on the basis of their plant communities.


Introduction
The analysis of landscape from the phytosociological point of view is a valuable tool for its comprehensive study, including its dynamism and heritage value. The common methodologies used in phytosociology are considered as an optimal choice in environmental management assessments of habitats, as has been recognized for decades [1][2][3][4][5][6][7][8]. It is based on floristic inventories of homogeneous areas and the evaluation of the taxa present according to their abundance and dominance. This method has proven very useful in obtaining knowledge on vegetation and its dynamics over increasingly large territories [9,10]. Despite the implicit subjectivity of such information, the enormous amount that has accumulated over the past century is currently viewed as an extraordinary database susceptible to statistical and multivariate analysis, using the inventory as a working unit. Researches have shown that these observations can be treated with a high degree of confidence [5,11,12]. The use of taxa and plant communities as indicators in land-use planning and their application in natural environment conservation policies is accepted in several countries, insofar as they are in themselves the object of such protection [13,14]. Some studies have applied the information contained in the study of vegetation (habitats) and their cartographic representation to territorial biological assessment criteria for natural areas [15][16][17]. Loidi (2008) [18] estimated the environmental or naturalistic value of natural

Rarefaction Index (RRF)
The results for the RRF index (Table S2) rendered four groups (ranges), considering the mean and the standard deviation ( Figure 2). The first range contained AMC and DEG, which had a low RRF and few communities; the second range contained TTR, DUA and EPU, which had a moderate RRF and also consisted of few communities; and the third range contained EST, which had a high RRF but few communities, and BMB, MIC, PES, CGN, MOD and MRF, which also had a high RRF as well as many communities. Finally, the fourth range contained BCA and DOÑ , which had a very high RRF and many communities.

Rarefaction Index (RRF)
The results for the RRF index (Table S1) rendered four groups (ranges), considering the mean and the standard deviation ( Figure 2). The first range contained AMC and DEG, which had a low RRF and few communities; the second range contained TTR, DUA and EPU, which had a moderate RRF and also consisted of few communities; and the third range contained EST, which had a high RRF but few communities, and BMB, MIC, PES, CGN, MOD and MRF, which also had a high RRF as well as many communities. Finally, the fourth range contained BCA and DOÑ, which had a very high RRF and many communities.

Priority Conservation Areas (PCA)
In the distribution (straight line) of PAs by Priority Conservation values (PCA; Table S1) and number of communities (Figure 3)

Priority Conservation Areas (PCA)
In the distribution (straight line) of PAs by Priority Conservation values (PCA; Table  S2) and number of communities (Figure 3), four groups of PAs were distinguished: DOÑ and BCA, with the highest PCA values and many plant communities; CGN and PES, which had relatively high values; MIC, MRF, EST, DUA, and AMC, which had moderate PCA values and small surface areas. BMB, DEG, EPU, MOD, and TTR presented PCA values of 0 because their plant communities are represented in other PAs; i.e., the 16 BMB and eight TTR communities were also represented in BCA (23); the 10 EPU and 21 MOD communities are also present in DOÑ (23); DEG is a very small area, with only six communities, which are also present in many of the other PAs.      Principal component analysis of the STD, RRF and PCA showed that component 1 has a variance of 91.8% and that rarefaction is the variable with the highest weight ( Figure 4 and Table 1). Five groups of PAs were clearly distinguished. Group A included DOÑ, BCA (closely linked by its communities to DOÑ) and CGN, in addition to PES (closely linked by its communities to CGN). Group B included EST and BMB, while Group C included MFR, MIC, and MOD, all of which had lower RRFs and STD, as they are closely linked to DOÑ. Group D included TTR and DUA, in addition to EPU. Finally, group E included AMC and DEG, both with a low number of communities (five and six, respectively), although in the case of AMC those present are almost exclusive to this area, hence its high STD and low RRF. In the case of DEG, the few communities present are found almost throughout the entire study area.

Sabinar.
Principal component analysis of the STD, RRF and PCA showed that component 1 has a variance of 91.8% and that rarefaction is the variable with the highest weight ( Figure  4 and Table 1). Five groups of PAs were clearly distinguished. Group A included DOÑ, BCA (closely linked by its communities to DOÑ) and CGN, in addition to PES (closely linked by its communities to CGN). Group B included EST and BMB, while Group C included MFR, MIC, and MOD, all of which had lower RRFs and STD, as they are closely linked to DOÑ. Group D included TTR and DUA, in addition to EPU. Finally, group E included AMC and DEG, both with a low number of communities (five and six, respectively), although in the case of AMC those present are almost exclusive to this area, hence its high STD and low RRF. In the case of DEG, the few communities present are found almost throughout the entire study area.

Correlation of the Legal Protection Index and DRA
The integration of the three indices into the new index DRA (acronym of STD, RRF and PCA) is shown in Table S1. We found a high correlation between the DRA and Legal Protection Index (LPI), which refers to the current protection status of a protected area (Section S1 and Table S1). However, it is worth noting the higher weight of the PCA and RRF with respect to it, with the STD being undervalued with the lowest correlation value ( Table 2).

Correlation of the Legal Protection Index and DRA
The integration of the three indices into the new index DRA (acronym of STD, RRF and PCA) is shown in Table S1. We found a high correlation between the DRA and Legal Protection Index (LPI), which refers to the current protection status of a protected area (Section S1 and Table S1). However, it is worth noting the higher weight of the PCA and RRF with respect to it, with the STD being undervalued with the lowest correlation value ( Table 2). The LPI and DRA values for each AP were distributed along a second-order polynomial trend line with R2 = 0.82 ( Figure 5). The values for LPI, STD, PCA and DRA are shown in Table S1.   The residuals obtained for each PA after fitting the curve obtained by comparing the LPI with the DRA allowed three categories to be distinguished, based on their absolute values and their graphical interpretation ( Figure 6): (1) those residuals whose absolute value is >1 (red bars) corresponded to PAs with a significant divergence between LPI and DRA (DEG, BCA, BMB and EPU); When clustering (Figure 7), the Phenon Line drawn at 80% similarity configured four groups, identified on the basis of DRA values, and consequently assignable to IUCN categories and/or those declared by the regional government of Andalusia (Spain) through the Andalusian Network of Natural Spaces [25].  When clustering (Figure 7), the Phenon Line drawn at 80% similarity configured four groups, identified on the basis of DRA values, and consequently assignable to IUCN categories and/or those declared by the regional government of Andalusia (Spain) through the Andalusian Network of Natural Spaces [25].

Discussion
The World Commission on Protected Areas (IUCN-WCPA) [26], defined "protected area" as: A clearly defined geographical space, dedicated and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values, while specifying among other principles that the fundamental objective should be the conservation of nature, as well as maintaining or even increasing the degree of naturalness of the ecosystem being protected. It also establishes seven categories in a non-hierarchical system, to be applied in the context of national protected area systems and as part of the ecosystem approach: Ia-Strict nature reserve; Ib-Wilderness area; II-National park; III-Natural monument or feature; IV Habitat/species management area; V-Protected landscape; VI-Protected areas with sustainable use of natural resources. Each of these categories was described on the basis of a number of criteria (distinctive values, role in the landscape or uniqueness), although it acknowledges that the assignment to any given category depends more on the intended use by the management authority rather than on any fixed and unchangeable set of criteria. Such categorization is important because it works as the basis of the definition of the management objectives. This is even more the casewhen the allocation of categories has traditionally been the responsibility of governments, which has led to this situation being questioned. In any case, the IUCN's own Guidelines for the Application of Protected Area Management Categories [26] recommend that "The only principle that should be applied in assigning categories is the appropriateness of the management objective assigned to the protected area within the system in relation to the ecological needs and threats to the species or ecosystem, in the context of the whole territory or marine environment in which the biodiversity occurs".
More recently, the criteria for determining the areas to be protected and their categorization evolved and advanced with the help of conservation biology, technological tools such as geographic information systems and statistical processing based on biodiversity databases. The quality and beauty of the landscapes are no longer the only criteria for the selection of an area; the representativeness and complementarity that a reserve offers for the protection of biodiversity are now incorporated [27].
The method of inventorying the plant communities should integrate broad ecosystemic information, and it should allow for hierarchization using the vegetation association as a basic unit. Such a hierarchy would include the association in successive higher levels of alliances, orders and classes that would allow a dynamic definition of the territory. For this, we propose an integrated index (DRA). This cumulative index strengthens three variables for the proposal of PA management: (1) Distinctness, understood as a measure of the diversity of a mosaic of ecosystems that distinguish a PA. It is increasingly recognized that appropriate measures of biodiversity within a particular taxonomic group should not merely be functions of the number of species and their relative abundance, but should also include information on the taxonomic relatedness of the species in question, which has led to the development of taxonomic distinctiveness (TD) [28]. On the same biological basis, syntaxonomy provides, through the hierarchization of communities and their treatment by the same method (STD), complementary information on the relationships of habitat diversity.
(2) Rarefaction, understood as the degree of rarity considering the extent of the PA. One of the challenges in determining the category of protection to be assigned to a territory is to resolve the size dependence and compare richness when the number of plant communities is not equal. To achieve this, one of the possible options [29] is the interpolation of species richness using the rarefaction technique, based on the shape of the species-abundance curve [30,31]. This traditional rarefaction was used to calculate the expected number of plant communities by reducing the samples to a standard size, i.e., interpolating into the same number of communities those areas with the lowest abundance [22].
(3) Prioritization, or understanding that the richness of communities confers a selective value in the process of choosing categories.
It should be noted that according to the results obtained in the clustering, in more than 70% of the cases there is an exact allocation of the categories with the values obtained from the DRA. Only four units differ in the value of the assigned category: PES is a Protected Site declared by the Regional Government, but its geographical proximity and degree of similarity in the cluster imply that it should be integrated as a Natural Park together with CGN, from which it is separated by the city of Almeria. Both BMB and EST are declared as Natural Parks; however, their coastal communities make them more similar to Natural Reserves. Finally, the low STD of EPU led to its integration as a Natural Monument, rather than a Natural Reserve as it appears in the RENPA.
The results of this study confirm that the use of the phytosociological method to describe plant communities in territories provides valuable information for the management of areas and their protection. Thus, statistical analysis using complementary diversity indices-such as STD, RRF and ACP-can be used to compare the status granted to PAs in territorial planning, using a new integrative index (DRA).

Materials
We selected the same 14 protected areas chosen by Salvo Tierra et al. (2020) [24] ( Figure 8). These areas are included in the Andalusian Network of Natural Spaces. In line with the findings of Pereña (2018) [25], we created a syntaxonomic scheme of the set of plant communities observed that included their corresponding alliance, order, and class (Appendix A). Based on this information, we constructed a Basic Data Matrix (BDM) of the plant communities in each protected area (Table S2), and assigned an identifying abbreviation to each syntaxon.
We selected the same 14 protected areas chosen by Salvo Tierra et al. (2020) [24] ( ure 8). These areas are included in the Andalusian Network of Natural Spaces. In line the findings of Pereña (2018) [25], we created a syntaxonomic scheme of the set of p communities observed that included their corresponding alliance, order, and class ( pendix A). Based on this information, we constructed a Basic Data Matrix (BDM) of plant communities in each protected area (Table S1), and assigned an identifying ab viation to each syntaxon.  Based on the results of the foregoing, numerical values were obtained for each Protected Area (Table 3).

Method
We used three statistical diversity indices to determine whether certain areas have been granted their correct level of protection, based on their plant communities: syntaxonomic distinctiveness (STD), rarefaction (RRF) and level of prioritization of conservation areas (PCA). These were multiplied in order to obtain the DRA index. This was applied to the protected areas on the Andalusian coast and compared with their current category of protection.
So far, these three indices have always been applied to conventional taxonomy, but in this case, we further applied them to syntaxonomy because just as a taxon carries a pool of ecological, biogeographical, evolutionary, etc. information, syntaxons provide equally significant information. In this case, this is relevant because they are communities of plant species that converge in the same territory with identical ecological conditions.

Syntaxonomic Distinctness (STD)
Taxonomic structures are based on the number of taxa included at the highest levels of the taxonomic hierarchy therefore, in the case of using species as a unit, the higher hierarchical levels would be genus, family, order, subclass/class. They have been widely used in ecology to measure impacts on ecosystem diversity [20,28,[32][33][34][35]. Among the estimators of taxonomic richness, taxonomic distinctness (TD [35] and references therein) is particularly significant because it is less influenced by sample size than by species diversity. Moreover, TD may be a more sensitive univariate index of community disturbance than species diversity [20]. Syntaxonomic structures provide relevant information, and thus have the same analytical potential as TD [36] for measuring impacts on territorial diversity. Clarke and Warwick (1998) [28] proposed the ∆+ index (Delta-plus) to estimate TD. We extend this index to the syntaxonomic distinction (STD) of PAs. This index measures the overall average length of the syntaxonomic path between two randomly chosen PAs.
The ∆+ index for STD would be calculated according to the following equation: where s is the number of phytosociological associations observed in a given PA, and w ij is the weight or distinction value given to each syntaxonomic branch of the hierarchical classification from association i to the first common node j (alliance), and so on up to the third and fourth level (order and class).

Rarefaction (RRF)
The rarefaction (RRF) index was used to compare observed diversity between locations that have not been sampled equally [16]. To make a fair comparison of these sites, subsamples (i.e., communities) were drawn from the larger sample, and the expected richness in each unit (i.e., PAs) was calculated based on abundance distributions in the larger sample. The process was repeated for subsamples of different sizes. The final value of the RRF index showed the expected value of community richness according to the sample size of each PA [21]. This correction allows for the direct comparison of the richness of two samples that were originally of different sizes.
Hurlbert (1971) [31] proposed the following algorithm to determine the RRF: where qi = ( n−zi n ) ( N n ) would represent the probability that community i does not appear in a sample of size n, z i is the number of i communities, and (N/n) is the binomial coefficient or the number of ways in which n out of N can be chosen. This algorithm was calculated by applying an analytical solution known as "Mao's tau", in which standard deviations are converted to 95% confidence intervals [37].

Priority Conservation Areas (PCA)
The method proposed by Vane-Wright et al. (1993) [38] was used to determine the value of the PCAs. In a first cycle, the PAs with the most communities were selected. Subsequently, the communities within the PA were eliminated from the BDM. This procedure was repeated in successive cycles that included the remaining communities that had not been included in the PAs already selected. The values obtained in each cycle refer to the total number of communities in the study area.
With the values of these three variables (STD, RRF and ACP), a Principal Component Analysis (using the variance-covariance matrix method) was implemented to check the influence of each one of them on the different protected areas, based on their vegetation communities, as well as the one that presented a higher level of incidence in the way these areas are grouped.
The calculations of these three indices were carried out using PAST v. 4.11 software [39,40].

Calculation of the DRA Index
The proposed DRA index was established as a tool for assessing the suitability of the assignment of protection categories to natural areas. For this reason, it was calculated by means of the product of the values obtained and subsequently standardized from the DST, RRF and PCA indices, all of which were independent. In this way, the range of variation of the values obtained was adjusted to a similar range to that obtained for the LPI (Table S1).
The values obtained from the calculation of the three indices and DRA are correlated (Spearman's non-parametric correlation) with the LPI, which is the reference index for the current protection status.
DRA values were compared with those of the LPI (Table S1) to check how they are distributed along the trend lines of both. This was done with excel V. 2204. With the residues obtained, divergences between the DRA and the LPI in the different protected areas were checked using the "confidence tunnel" method [28].
Finally, in order to test the functionality of the new DRA index, a clustering of the PAs studied was carried out on the basis of the three components of the index and the current protection categories (LPI), taking into account the IUCN categories for the protection of areas and that of the Andalusian Network of Protected Spaces, and drawing a Fenon Line at 80% similarity. This was calculated with the PAST 4.1.1 software using the Euclidean distance and the WPGMA algorithm (co-phenetic correlation of 0.87).

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