*Article* **Ecological and Biological Properties of** *Satureja cuneifolia* **Ten. and** *Thymus spinulosus* **Ten.: Two Wild Officinal Species of Conservation Concern in Apulia (Italy). A Preliminary Survey**

**Enrico V. Perrino 1,\*, Francesca Valerio 2, Shaima Jallali 1, Antonio Trani <sup>1</sup> and Giuseppe N. Mezzapesa <sup>1</sup>**


**Abstract:** This study evaluated the effects of ecology (plant community, topography and pedology), as well as of climate, on the composition of essential oils (EOs) from two officinal wild plant species (Lamiales) from Apulia, namely *Satureja cuneifolia* Ten. and *Thymus spinulosus* Ten. Few scientific data on their chemical composition are available, due to the fact that the first has a limited distribution range and the second is endemic of southern Italy. Results for both species, never officially used in traditional medicine and/or as spices, showed that the ecological context (from a phytosociological and ecological point of view) may influence their EO composition, and hence, yield chemotypes different from those reported in the literature. *S. cuneifolia* and *Th. spinulosus* can be considered good sources of phytochemicals as natural agents in organic agriculture due to the presence of *thymol* and *α-pinene*. Overall, the obtained trend for EOs suggests a potential use of both species as food, pharmacy, cosmetics and perfumery. Hence, their cultivation and use represent a positive step to reduce the use of synthetic chemicals and to meet the increasing demand for natural and healthier products.

**Keywords:** correlation; ecology; essential oils; Lamiaceae; vegetation

#### **1. Introduction**

Medicinal and aromatic plants (MAPs) have a long history of being used for multiple purposes, including for food and therapeutics. Before the advances of modern medicine, ancient civilizations discovered a wealth of useful therapeutic agents from within the plant and fungi kingdoms [1]. Knowledge of these medicinal preparations was passed down through generations and was occasionally recorded in the herbal literature [2]. Hence, the use of these wild officinal plants became a significant aspect of populations' cultural heritage and transformed into different traditions that kept on being passed down [3,4]. Much research in the field of biology, chemistry and medicine is directed at the identification and characterization of plant secondary metabolites with a pharmacological activity, which may be candidates for the synthesis of new drugs [5].

MAPs are defined as plants possessing aromatic and/or medicinal characteristics that can be extracted in the form of essential oils (EOs) through their secondary metabolites of bioactive properties [6]. These extracts are synthesized for several reasons, including plant protection against pathogens, insects, pollinators attraction and for allelopathic activity [7–9].

EOs are natural complexes of volatile and semi-volatile organic compounds, responsible for specific aroma, flavor and fragrance of certain plant species [10]. They can be stored by all plant organs (flowers, leaves, stems, twigs, seeds, fruits, roots, wood, or bark) and are synthesized in different histological structures such as secretory cells, cavities, canals or glandular trichomes [8,11].

**Citation:** Perrino, E.V.; Valerio, F.; Jallali, S.; Trani, A.; Mezzapesa, G.N. Ecological and Biological Properties of *Satureja cuneifolia* Ten. and *Thymus spinulosus* Ten.: Two Wild Officinal Species of Conservation Concern in Apulia (Italy). A Preliminary Survey. *Plants* **2021**, *10*, 1952. https:// doi.org/10.3390/plants10091952

Academic Editor: Calvin O. Qualset

Received: 25 August 2021 Accepted: 14 September 2021 Published: 18 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

5

The chemical structure of plant EOs falls into two distinct classes: terpenes and phenylpropanoids. Terpenes are formed by the combination of isoprene units such as monoterpenes (two units), sesquiterpene (three units), diterpenes (four units), and phenylpropanoids, which are aromatic and non-terpenoid compounds responsible in part for the fragrance of the plant [11,12]. These volatile oils represent a wide chemical diversity and are responsible for many plant applications, such as antimicrobial, antiparasitic, insecticidal, herbicidal, antioxidant and other medicinal properties.

The Lamiaceae family (formerly known as Labiatae) is the largest family of the order Lamiales, containing more than 245 genera and 7886 species [13], which grow in different natural ecosystems worldwide, with a particularly high concentration in the Mediterranean region, including Italy, providing a great level of biodiversity [14]. Most of the species belonging to this family are aromatic and possess EOs, making them valuable in cosmetics, perfumery, agriculture and medicine [15,16].

EOs can present a wide variability in their yield and chemical composition due to a series of intrinsic and extrinsic factors that affect the plants in their environment. In fact, it is well documented that the same species, under different environmental and geographical conditions, can produce EOs with different chemical profiles and biological properties [17–19]. In addition, some of these taxa are rare or endemic and have not yet been investigated for their ecology, chemical composition and biological properties.

The aim of the present research was to investigate two wild species of conservation concern (SCC; species for which we have concerns about its ability to remain on a landscape for a long time), *Satureja cuneifolia* Ten. (wild savory) and *Thymus spinulosus* Ten. (little thorn thyme), in order to: (a) perform ecological studies to understand the relationships existing between plant associations (sociology) and their surrounding habitat; (b) verify and eventually explain the ecological context influences on the composition of EOs; (c) determine chemical composition; (d) evaluate the potential for commercial purpose.

#### **2. Study Area**

Samplings were carried out in 2020 in Apulia. For *Satureja cuneifolia*, two sites, San Basilio (Mottola, province of Taranto) and Difesa di Malta (Fasano, province of Bari), were detected, while for *Thymus spinulosus*, three sites were identified, one of which, San Basilio (Mottola), was the same as *S. cuneifolia*, and the other two sites being Scannapecora (Altamura, province of Bari) and S. Egidio in S. Giovanni Rotondo (province of Foggia) (Figure 1).

The investigated sites have a Mediterranean macro-bioclimate, except S. Egidio, which is cooler, resulting Temperate. The bioclimate is oceanic, particularly the semicontinental transition at S. Egidio and Scannapecora, while the ombrotype is always the subhumid and dry type at Difesa di Malta (Table 1) (Phytoclimatic map—Italian Ministry for the Environment, Land and Sea, http://www.pcn.minambiente.it/viewer, accessed on 26 May 2021) [20] (Table 1).


**Table 1.** Phytoclimate at the investigated sites of *S. cuneifolia* (*Sc*) and *Th. spinulosus* (*Ts*).

**Figure 1.** Site locations of *Satureja cuneifolia* and *Thymus spinulosus*.

Within the San Basilio site in the province of Taranto, the most widespread pedotype is Lithic Ruptic–Inceptic Haploxeralf fine, predominantly clayey, very thin, and very rocky with substrate within 50 cm [21]. At Scannapecora in the province of Bari, the geological type is that of Skeletal limestones of neritic and carbonate platform facies (Upper Cretaceous), on Tyrrhenian carbonate reliefs with material defined by calcareous sedimentary rocks and climate from Oceanic to Suboceanic Mediterranean, partially mountainous. The Difesa di Malta and S. Egidio sites share the same ecopedology as the Scannapecora site but different geopedology and geolithology, respectively, with Terrigenous-skeletal limestones such as "Panchina" (Pleistocene) and Colluvial, terraced alluvial, fluviolacustrine and fluvioglacial deposits (Pleistocene) (Geological, Geolithological and Ecopedologic maps—Italian Ministry for the Environment, Land and Sea, http://www.pcn.minambiente.it/viewer, accessed on 26 May 2021) [20] (Table 2).

**Table 2.** Geopedology, geolithology and ecopedology in the investigated sites of *S. cuneifolia* (*Sc*) and *Th. spinulosus* (*Ts*).


#### **3. Materials and Methods**

#### *3.1. Vegetation Analysis*

The field inspections of *Satureja cuneifolia* and *Thymus spinulosus* were carried out in 2020. A total of five vegetation surveys were conducted for both species, following the phytosociological method of the Zurich–Montpellier Sigmatist School [22], updated by Géhu and Rivas-Martínez [23], using the coefficient of abundance–dominance of identified species, well known to field botanists, and concerning the covering percentage of the relief area (value 5: species with more 75% coverage; value 4: species with coverage from 50% to 75%; value 3: species with coverage from 25% to 50%, value 2: species with coverage from 5% to 25%: value 1: species with coverage from 1% to 5%; value +: species with coverage less than 1%). Identification code, taxon, location, date, and geographical position (expressed in WGS84—World Geodetic System 1984) are reported in Table 3, while phytosociological and topographical data (identification code, altitude (m a. s.), aspect, slope (◦), relevé area (m2), stoniness (%), rockiness %), cover total (%), soil depth detected (cm), macroclimate, geolithological and eco-pedological types), endemicity, number of species identified, and number of individuals collected for laboratory analyses are reported in Tables 4 and 5. Other identified and recorded species not used for the phytosociological classification are not reported here. A total of 5 specimens for each species and for each reléve were collected and deposited at the official herbarium of Bari University (Italy) (*Herbarium Horti Botanici Barensis*—BI).

**Table 3.** Location, date retrieved, and geographical position in the investigated sites of *S. cuneifolia* (*Sc*) and *Th. spinulosus* (*Ts*).


**Table 4.** *Phagnalo saxatilii-Saturejetum cuneifoliae* Biondi & Guerra 2008.


**Table 4.** *Cont.*


*Identification Code* **Sc1 Sc2** ††—Tyrrhenian carbonate reliefs with material defined by calcareous sedimentary rocks and climate from Oceanic to Suboceanic Mediterranean, partially mountainous **Macroclimate** =—Mediterranean ==—Temperate **Bioclimate** ♦—Mediterranean Oceanic ♦♦—Oceanic ♦♦♦—Oceanic-semicontinental transition **Ombrotype** #—Dry ##—Subhumid

**Table 4.** *Cont.*

Identification of the taxa was carried out according to Flora d'Italia [24] and Flora Europea [25], with nomenclature standardized by "An updated checklist of the vascular flora native to Italy" [26] and "An updated checklist of the vascular flora alien to Italy" [27], while the syntaxonomic framework was conceived by several contributions [28–30], reported in the phytosociological tables and summarized in the syntaxonomic scheme.

The aerial parts of the studied plant species were harvested in the same sites and date of vegetation surveys. The number of individuals collected (from 15 to 100) for the laboratory analysis is reported in the phytosociological tables (Tables 4 and 5).

**Table 5.** Physico-chemical characteristics of *S. cuneifolia* (*Sc*) and *Th. Spinulosus* (*Ts*) for each investigated soil site.


#### *3.2. Essential Oils Analysis*

The essential oils (EOs) of the two wild plant species were extracted by hydrodistillation [31] using a Clevenger type apparatus for 4–5 h. Air-dried plant material (100 g) was covered with 500 mL of distilled water in a 1L volume distillation flask (extraction ratio

10

1:5 w/vol). EOs were collected in amber glass vials, weighed, and stored at 4 ◦C. EO yields (% *w*/*w*) were determined as grams of EOs per 100 g of dry weight plant material.

Identification of the compounds present in the EO extracts was carried out by gas chromatography coupled with mass spectrometry (GC–MS) using a Clarus 680 GC interfaced with a Clarus SQ8C single quadrupole mass spectrometer (PerkinElmer, Waltham, MA, USA) equipped with an Elite-5 MS (PerkinElmer, Waltham, MA, USA) fused silica capillary column (30 m × 0.25 mm and 0.25 μm film thickness). Mass spectra of the target compounds were obtained by an electron impact ionization system with standardized ionization energy of 70 e V. Helium 5.5 was used as carrier gas at a constant flow rate of 1ml/min. Mass transfer line and injector temperatures were set at 280 ◦C and the oven temperature was programmed from 50 ◦C to 160 ◦C at 5 ◦C/min, then raised to 250 ◦C at 10◦ C/min, and held at the final temperature for 5 min. Diluted samples (1:10, *v*/*v*, in hexane) were injected in split mode with a split ration of 1:100. Data were collected in full scan mode in the range 40−300 amu. A solvent delay of 4 min was applied. Qualitative results include compound identification and area percentage of related peaks in the total ions chromatogram (Table 5).

Compounds identification was performed by both Kovats retention indexes (RI) [32,33] and mass spectra (MS) searches in NIST and Wiley databases (Table 6).

**Table 6.** EO composition as area percentage of *S. cuneifolia* (*Sc*) and *Th. spinulosus* (*Ts*) with the identification method used. RT: retention time; RI: Kovats retention index; RI: experimental retention index; RI ref.: bibliographic retention index according to [32].



**Table 6.** *Cont.*

The chemical composition (%) of the main compounds of *Satureja cuneifolia* and *Thymus spinulosus* is reported in Table 7. Statistical analysis was not applicable due to the low number of samples per species and site, limited amount of material per sample, and few and variable numbers of individuals per specimen (15–70), also due in part to the fact that the study was carried out on two species of conservation concern, which limits the harvest of material.

**Table 7.** *Hippocrepido glaucae-Stipion austroitalicae* Forte & Terzi in Forte, Perrino and Terzi 2005.


**Table 7.** *Cont.*


**Table 7.** *Cont. Identification Code* **Ts1 Ts2 Ts3** ††—Tyrrhenian carbonate reliefs with material defined by calcareous sedimentary rocks and climate from Oceanic to Suboceanic Mediterranean, partially mountainous **Macroclimate** =—Mediterranean ==—Temperate **Bioclimate** ♦—Mediterranean Oceanic ♦♦—Oceanic ♦♦♦—Oceanic-semicontinental transition **Ombrotype** #—Dry ##—Subhumid

#### *3.3. Soil Analysis*

Soil samples were collected from each vegetation sample of topsoil in the following way: in 5–9, subsamples (depending on the habitat extension; each subsample ∼ 100 g) were taken from the first 0–20 cm of depth.

Soil samples were sieved to <2 mm and stored for physical and chemical analysis. Particle size distribution was determined by the pipette method, whereas textural class was categorized according to the USDA classification [34]. Soil pH was measured both in water and saline solution using 1:2.5 ratio (*w/v*). Electrical conductivity (EC) was measured on a soil to distilled water ratio (1:2, *w/v*). Soil organic carbon (OC) was determined according to the Walkley and Black method as described by Nelson and Sommers [35]. Total nitrogen (N) was measured following the Kjeldahl method as described in Bremner [36]. Available P was measured in sodium bicarbonate alkaline soil extracts and determined colorimetrically [37]. The total carbonate content in the soils was determined using a Dietrich–Fruhling calcimeter. Exchangeable bases were extracted by BaCl2 (1 M) and triethanolamine solution, while available microelements were extracted by DTPA and triethanolamine solution. After extractions, both were analyzed by ICP-OES.

#### **4. Results and Discussions**

#### *4.1. Satureja Cuneifolia Ten.*

*S. cuneifolia* (synonym: *Satureja montana* subsp. *cuneifolia* (Ten.) O. Bolós & Vigo) is shown in Figure 2.

*Herbarium* samples are shown in Figure 3: The sample from San Basilio, in the municipality of Mottola (Taranto), garrigue, 26 May 2020, *legit* and *determinavit* determinative *E.V. Perrino* 42457 (*Herbarium Horti Botanici Barensis*—BI), is shown in Figure 3a. The sample from Difesa di Malta, in the municipality of Fasano (Brindisi), garrigue, 3 June 2020, *legit* and *determinavit E.V. Perrino* 42458 (*Herbarium Horti Botanici Barensis*—BI), is shown in Figure 3b.

**Figure 2.** *S. cuneifolia* (**a**) in flowering; (**b**) in its habitat, *Phagnalo saxatilii-Saturejetum cuneifoliae*. San Basilio, 26 March 2020.

**Figure 3.** *S. cuneifolia* herbarium samples: (**a**) San Basilio (Mottola) (BI 42457); (**b**) Difesa di Malta (Fasano) (BI 42458).

The genus *Satureja* L. comes from the Latin "*satureia*" and was first named by Roman writer Pliny [38]. It means "*herb of satyrs*" and, for this reason, its cultivation was banned in monasteries [39]. *S.* belongs to the Lamiaceae family, sub-family *Nepetoideae*, and the tribe *Mentheae* [40] including about 200 wild species of herbs and shrubs, often aromatic, widely distributed in the Mediterranean area, Asia and boreal America [41,42]. More than 30 of them grow in eastern parts of the Mediterranean area [43,44], up to 1200 m above sea level, in arid, sunny, stony and rocky habitats [45], 5 of which are in Italy: *S. cuneifolia* Ten., *S. montana* L. subsp. *montana*, *S. montana* L. subsp. *variegata* (Host) P.W. Ball, *S. subspicata* Bartl. ex Vis. subsp. *liburnica* Šilic, and *S. thymbra* L. [26]. *S. cuneifolia* Ten. is a SE European entity reported in Italy, Albania, Greece, Croatia and Montenegro [46,47], with its western limit of distribution in Italy, where the species occurs only in Apulia, throughout the Region, in Salento, Murge and Gargano [48–52], and in the nearby Basilicata, in the Murge di Matera [53]; it was recorded by mistake at Laino Borgo in Calabria [24,26]. As highlighted by some authors [54], this taxon is one of the many plant species with eastern distribution that exhibit a disjunct Italian distribution more or less restricted to the Apulian region, such as *Bromus parvispiculatus* H. Scholz [55], *Carex phyllostachys* C. A. Mey. [56,57], *Cerinthe retorta* Sm. [58], *Linum elegans* Spruner ex Boiss. [59], *Scrophularia lucida* L. [60,61], and *Ophrys oestrifera* complex [62].

*S. cuneifolia* grows from sea level to 600 m of altitude, and is a diagnostic taxon of three plant communities: (a) *Phagnalo saxatilii-Saturejetum cuneifoliae* Biondi & Guerra (2008) associated with *Phagnalon saxatile* (L.) Cass.; (b) *Sedo ochroleuci-Saturejetum cuneifoliae* Di Pietro & Misano (2010) associated with the endemism *Petrosedum ochroleucum* (Chaix) Niederle subsp. *mediterraneum* (L. Gallo) Niederle; (c) *Stipo austroitalicae-Seslerietum juncifoliae* Di Pietro & Wagensommer (2014) subass. *typicum* [54]. These vegetations fall into the *Ononido-Rosmarinetea* class, with the first two (a, b) described for central of Apulia in the "Gravine" gorges [49,63], while the third (c) for the Gargano [50,54]. The *Stipo austroitalicae-Seslerietum juncifoliae* subass. *typicum* grows on small rocky outcrops emerge from a matrix composed of limestone debris, found within the steep slopes of the gorges, which cut across the southern side of the Gargano promontory [54]; it could be considered a distinct aspect of the habitat "Calcareous rocky slopes with chasmophytic vegetation" (code 8210) according to Directive 92/43 EEC [64,65].

*Satureja* species have been traditionally used in the treatment of many diseases such as nausea, indigestion, cramps, diarrhea, infectious diseases and muscle pains [66,67], in relation to the EOs secreted by the glands found on their leaf surface, as with other aromatic plants belonging to some genera of the Lamiaceae family.

Previous studies from the Mediterranean area showed variability in chemical composition in EOs of *S. cuneifolia*, attributed to many factors such as climate, regional, local and ecological conditions, growth stages, harvesting period, and genetics [67–71]. In three different localities of Dalmatia (Croatia), corresponding to as many stages of development, the yield of the EOs ranged from 0.1% to 0.6%, and all samples had a low percentage of thymol (0–3.5%) and carvacrol (0–16.3%), but were relatively rich in α-pinene (1.3–20.7%), limonene (1.8–17.4%), γ-terpinene (0–5.6%), and p-cymene (1.8–14.8%) [69]. In a research study carried out in the southwestern part of Montenegro (Lov´cen National Park), a high percentage of linalool (20.3%) has been found, the most abundant compound, followed by much lower values of α-terpineol (3.8%) and borneol (3.6%), while limonene (1.1%) and α-pinene (0.7%) were under-represented, and carvacrol and thymol were even absent [67]. In high altitude conditions and precisely from Sogut mountain in Turkey, *S. cuneifolia* collected during the flowering stage in August recorded elevated values of carvacrol (45%), p-cymene (21.6%), and then, lower but appreciable values thymol (9%), γ-terpinene (4.3%), and borneol (2.5%) [72]. A further study [70] conducted in 12 wild sites far east of the Italian peninsula at lower altitudes than Turkey showed significant percentage changes in the main compounds observed: α-pinene (0.3–26.4%), linalool (0.2–45.8%) and borneol (1.2–41.1%), with linalool and α-pinene not always present. An accurate reading of the Italian data also highlights that linalool and α-pinene reached high values (>35% and >14%, respectively) only in four samples, while borneol was the most present chemotype.

*S. cuneifolia* from Croatia and Montenegro is very different from that collected in the Mediterranean parts of Turkey and eastern Italy (Apulia). EOs of Croatian and Montenegro origin have a low content (even absent) of both thymol and carvacrol, as observed in Italian sites, while, in general, EOs of plants from Turkey have a mean percentage of thymol accompanied by a high percentage of carvacrol or vice versa. In Italy, the values are even more extreme and variable because not only have thymol and carvacrol not been recorded, but they have considerable fluctuations within the same compound in comparable environments. The authors [70] explain that this intraspecific polymorphism indicates that homogeneous environments may host different genotypes.

The data collected about *S. cuneifolia* show that the vegetation belongs to *Phagnalo saxatilii-Saturejetum cuneifoliae* at both investigated sites (Table 4), where there is the same Mediterranean macroclimate, but different bioclimate and ombrotype, respectively: Dry Mediterranean Oceanic under the influence of the sea at Difesa di Malta, and Subhumid Oceanic in the hinterland at San Basilio. In addition to the climatic differences, there are also geolithological and ecopedological differences, with a more recent Pleistocene geology near the sea, where we find Terrigenous-skeletal limestones such as "Panchina", and older in the hinterland (Upper Cretaceous) with skeletal limestones of neritic and carbonate platform facies. Differences were confirmed by Hilly reliefs with undifferentiated tertiary sedimentary rocks at Difesa di Malta and Tyrrhenian carbonate reliefs with material defined by calcareous sedimentary rocks at San Basilio. This means that the two sites only have the same macroclimate, but different bioclimate, geology, and pedology characteristics, and partially different plant communities. In fact, within the same plant association, we observe less human-made interference in San Basilio than Difesa di Malta, and a greater biodiversity due to the many species referable to the *Hippocrepido glaucae-Stipion austroitalicae* alliance and *Scorzoneretalia villosae* order, enhanced with two endemic taxa: *Scorzonera villosa* subsp. *columnae* and *Thymus spinulosus* Ten. The vegetation detected in San Basilio deserves more analysis from a syntaxonomic point of view, through ad hoc surveys in the surrounding areas, to better define this distinct aspect with *S. cuneifolia*.

The soil parameters show some differences between the two sites (Sc1, Sc2). In the more disturbed aspect, less natural, in a cultivated context, as Difesa di Malta (Sc1), a fine-silt loam soil has been recorded, characterized by poor total carbonate, very high phosphorus availability and high total nitrogen and organic carbon, while at San Basilio (Sc2), in a more sandy soil, lower and moderate values were recorded, except for the total carbonate which, here, has a much higher concentration (Table 5), confirming the geopedological and geolithological features (Table 4).

The exchangeable cation and available microelement (macro and micronutrients, respectively) soil contents showed that all samples are rich in exchangeable cation but less rich in available microelements, with differences in types and percentage at both sites. The available microelements have comparable concentrations at both sites, except for Zn and Mn, which are more abundant in natural contexts, unlike the Fe that prevails in agricultural land (Table 5).

The environmental, geological, ecological, vegetational and pedological differences have a slight effect on the phytochemical properties of *S. cuneifolia*, remembering that the San Basilio site is located at a higher altitude (274 m), has a major plant biodiversity, and older soils, compared to the maritime site of Difesa di Malta (49 m), which is situated in a more agricultural context and in a more recent soil. In particular, the correlation between EO compositions was meaningful and positive with the monoterpenoid alkene and monoterpenoid alkene alcohol, which are the most abundant classes in *S. cuneifolia* EOs. Additionally, monoterpenoid alkenes have a significant positive correlation with soil pH, and significant negative correlation with total nitrogen and soil organic carbon (Table 5).

In all, 36 compounds were identified at the Difesa di Malta (Sc1) and San Basilio (Sc2) sites, but with some quantitative differences (Table 6). Most of the compounds are present with similar percentages at both sites. In particular, α-Pinene, the most abundant phytochemical, is present with the highest percentage at both sites (36.8% and 38.8%, respectively), the same slight difference occurs for limonene (5.1% and 6.4%), but significant differences arise for α-Terpineol (11% and 17.1%), trans-3-caren-2-ol (4.1% and 1.9%), and Borneol (6.9% and 4.1%). The rest of the compounds for both sites did not exceed 5%.

Finally, even other minor compounds, not mentioned, were present with slightly different percentages.

#### *4.2. Thymus Spinulosus Ten.*

*Th. spinulosus* (synonym: *Th. paronychioides* Celak.) is shown in Figure 4.

**(a) (b) Figure 4.** *Th. spinulosus* (**a**) in flowering; (**b**) in its habitat, *Hippocrepido glaucae-Stipion austroitalicae*. Scannapecora, 03.06.2020.

*Herbarium* samples are shown in Figure 5: The samples from San Basilio, in the municipality of Mottola (Taranto) (garrigue, 5 June 2020, *legit* and *determinavit E.V. Perrino* 42454 (*Herbarium Horti Botanici Barensis*—BI) and San Egidio, in the municipality of San G. Rotondo (Foggia), garrigue, 12 June 2020, *legit* and *determinavit E.V. Perrino* 42455 (*Herbarium Horti Botanici Barensis*—BI) are shown in Figure 5a. The samples from Scannapecora, in the municipality of Altamura (Bari), garrigue, 3 June 2020, *legit* and *determinavit E.V. Perrino* 42456 (*Herbarium Horti Botanici Barensis*—BI) are shown in Figure 5b.

The genus *Thymus*, described by Carl Linnaeus in *Species Plantarum* [73], is one of the most important genera within the family Lamiaceae, due to its high number of species, commercial uses and medicinal features [74]. It includes 220 accepted taxa distributed in Europe, Northwest Africa, Ethiopia, Asia and southern Greenland [75–77]. The center of diversity of the genus is the Mediterranean region, where the typical endemism grows especially in the west sector (Iberian peninsula, Northwest Africa) [78]. In Italy, 20 species of *Th.* are reported, with 5 endemims of the southern regions (*Th. paronychioides* Celak., ˇ *Th. picentinus* (Lacaita) Bartolucci, *Th. praecox* Opiz subsp. *parvulus* (Lojac.) Bartolucci, Peruzzi & N.G. Passal., *Th. richardii* Pers. subsp. *nitidus* (Guss.) Jalas, and *Th. spinulosus* Ten.) [26], where they occur mainly in rocky habitats, from hills to mountain tops. *Th. spinulosus* is reported in Campania, Apulia, Basilicata, Calabria and Sicily, while its presence is uncertain in Molise, no longer recorded in Latium, and reported by mistake in Abruzzo [26]. It grows in arid stony slopes, breached, clearings of deciduous oaks, from sea level up to about 1100 m of altitude and flowering in May–June. In Apulia, it is widespread in the perennial grasslands, especially in the pseudo-steppes of *Stipa austroitalica* Martinovský [29,54,79], which is considered is habitat (Eastern sub-Mediterranean dry grasslands (*Scorzoneretalia villosae*)" (code 62A0)) under Directive 92/43 EEC [64,65], while it is very rare in typically rocky vegetation [54].

Many species of *Thymus* are popular in the traditional medicine of many countries and nations as a source of valuable crude drug, where they are used for acute and chronic bronchitis, and the leaf extracts are prescribed as expectorant and analgesic [80]. In medicine and perfumery, *Th. serpyllum* and *Th. vulgaris* are widely used [81,82], and in general, many species of this genus show antimicrobial, anti-inflammatory, antioxidant, cytotoxic, spasmolytic, and antinematodal activities. *Thymus linearis* Benth. is used in Kashmir regions for gastrointestinal problems and fever issues [83]. *Thymus capitatus* (L.) Hoffmanns & Link (=*Thymbra capitata* (L.) Cav.) were used in Southern Italy to make fumigations in the treatment of colds together with dried figs and chestnuts [84]. Traditionally, the plants

of the genus *Th.* are used as a source of EOs and are promising in pest management and control of harmful insects, especially thanks to thymol and carvacrol content, as well as citral, geraniol and nerolidol [85–87]. *Th.* still need to be explored as supplementary sources of botanical insecticides, to face the growing concerns arising from the use of chemical pesticides [88–90], with special reference to insecticide resistance [91]. *Th. spinulosus* is a less explored taxon for both phytochemistry and biological activities [92], and there are few works in the literature detailing information on its use in the kitchen such as the preparation called "sanguinaccio" in the Campania region, based on pig blood with dried fruit, sultanas and cooked wine [93]. The Sicilian data, obtained from four populations located above 1000 m of altitude, but in different soil conditions (siliceous and calcareous), gave oils with different composition, 62 components in total, which is inconsistent with the previous chemical composition of EOs described in the literature for this endemic taxon, suggesting a new chemotype, characterized by myrcene-limonene among monoterpenes and γ-muurolene, caryophyllene and germacra-1(10),4-dien-6-ol among sesquiterpenes as the main constituents [87].

**Figure 5.** *Th. spinulosus* herbarium samples: (**a**) San Egidio (S. G. Rotondo) (BI 42455); (**b**) Scannapecora (Altamura) (BI 42456).

The macroclimate in the two eastern sites (Ts1, Ts2) is Mediterranean, while in the Gargano, it is colder and temperate (Ts3). The bioclimate is Oceanic-semicontinental transition and subhumid Ombrotype in the two sunniest sites placed at higher altitudes (Ts2, Ts3), and Dry Mediterranean Oceanic type in cooler expositions, as observed at San Basilio (Ts1). The Gargano site differs from the other two sites also due to its geolithological characteristics, more recent (Pleistocene) for colluvial, terraced alluvial, fluviolacustrine and fluvioglacial deposits (Ts3) compared to those of the upper Cretaceous corresponding

to skeletal limestones of neritic and carbonate platform facies (Ts1, Ts2). Additionally, the geopedology is different, with Tyrrhenian carbonate reliefs and material defined by calcareous sedimentary rocks on Gargano (Ts3), and Hilly reliefs with undifferentiated tertiary sedimentary rocks in the other two sites (Ts1, Ts2). This means that the three sites only have the same macroclimate, but the Gargano site (Ts3) differs from the other two due to different bioclimate, geology, and pedology characteristics, and partially different plant communities.

The phytosociological data (Table 7) show within the same alliance *Hippocrepido glaucae-Stipion austroitalicae*, there are two different plant associations, *Acino suaveolentis-Stipetum austroitalicae* (Ts1, Ts2) and *Sideritido syriacae-Stipetum austroitalicae* (Ts3). At S. Egidio (Ts3), in a good environmental condition, in the absence of human disturbances, with 626 m altitude, 12◦ slope, a medium level of stoniness (10%) and a high level of rockiness (45%), and 85% coverage of 64 plant species, we observed within the *Sideritido syriacae-Stipetum austroitalicae* association identified by *Sideritis italica*, a richness of species of *Hippocrepido glaucae-Stipion austroitalicae*, *Scorzoneretalia villosae* and *Festuco-Brometea* communities, with a good coverage of *Stipa austroitalica* subsp. *austroitalica*, *Scorzonera villosa* subsp. *columnae*, both endemics, and *Phleum hirsutum* subsp. *ambiguum*, in addition to other species such as *Petrorhagia saxifraga* subsp. *gasparrinii*, *Hippocrepis glauca*, *Teucrium capitatum* subsp. *capitatum*, *Convolvulus cantabrica*, *Koeleria splendens*, and *Eryngium amethystinum*. Some transgressive species of *Helianthemetea guttati* Rivas Goday & Rivas-Mart. 1963 and *Lygeo sparti-Stipetea tenacissimae* Rivas-Mart. 1978 classes show natural catenal contact with these other plant communities. In the two easternmost sites (Ts1, Ts2), always in natural conditions, but at different topographic conditions, such as altitudes (Ts1 = 274 m, Ts2 = 561 m), exposition (N–NE, S–SW) slope degree (2◦ and 10◦), level of stoniness (5% and 20%) and rockiness (60% and 40%), we detected the same *Acino suaveolentis-Stipetum austroitalicae* association, characterized by three species, two of which are endemic and present at both sites (*Thymus spinulosus* and *Euphorbia nicaeensis* subsp. *japygica*), while *Clinopodium suaveolens* was present only at Scannapecora (Ts2). The taxa detected are comparable to those of Gargano (Ts3), except for the characteristics of the Festuco-Brometea class, which are less represented at both sites, explained by the lower plant biodiversity (Ts1 = 38 taxa, Ts2 = 44 taxa) than Gargano (Ts3 = 64 taxa), and only for Scannapecora (Ts2), also with good coverage of Ononido-Rosmarinetea transgressive species, which this plant association makes a mosaic.

The GC–MS results showed that EO composition consists of a total of 27 identified compounds for all three sites (Table 6). The environmental differences are reflected only partially in the chemical composition of EOs. Only low differences were observed on the abundance and patterns among the three sites, since they share 25 compounds; only 1 (Hotrienol) was exclusive to Ts1 (Hotrienol) and 1 to Ts2 (Caryophyllene oxide). The most interesting phytochemicals, having the highest presence in all samples, were thymol, p-Cymene and β-Ocimene. In particular, at San Basilio (Ts1), in the same climatic condition as Altamura (Ts2), but with lower altitudes and different topographic, geological and geopedological conditions, 28 compounds were detected, with the most abundant ones similar to the Altamura (Ts2) site. Slight differences in percentages were observed: thymol (42.9% and 48.8%, respectively), followed by p-Cymene (17.9% and 17.5%) and β-Ocimene (15.4% and 11.7%). The other compounds did not exceed 3% and 4%.

At the Gargano (Ts3) site, in a different climatic, environmental and partly different plant community, 28 compounds were extracted, comparable anywhere with the results of the other two sites (Ts1 and Ts2): thymol (45.9%), p-Cymene (17.5%) and β-Ocimene (11.0%). The abovementioned data differ from the Sicilian populations [86], and highlight a different chemotype, probably due to the lower altitude compared to the Sicilian stations, which exceed 1000 m above sea level, and maybe for the different latitude, while it seems that the climate, the pedological and lithological characteristics of the soil and the vegetation do not play a relevant role. In any case, the qualitative difference seems to be particularly relevant, since Sicilian populations show a number of components almost double that of Apulia.

#### *4.3. Syntaxonomical Scheme of Surveyed Vegetation*

*Festuco-Brometea* Br.-Bl. et Tx. ex Soó 1947

*Scorzoneretalia villosae* Kovaˇcevi´c 1959

*Hippocrepido glaucae-Stipion austroitalicae* Forte et Terzi 2005 in Forte, Perrino and Terzi 2005

*Hippocrepido glaucae-Stipienion austroitalicae* Biondi and Galdenzi 2012

*Acino suaveolentis-Stipetum austroitalicae* Forte et Terzi in Forte, Perrino et Terzi 2005 *Sideritido syriacae-Stipetum austroitalicae* Biondi and Guerra 2008

*Ononido-Rosmarinetea* Br.-Bl. in A. Bolòs y Vayreda 1950

*Cisto-Micromerietalia julianae* Oberd. 1954

*Cisto cretici-Ericion manipuliflorae* Horvati´c 1958

*Phagnalo saxatilii-Saturejetum cuneifoliae* Biondi and Guerra 2008

#### **5. Conclusions**

The composition of *Satureaja cuneifolia* EOs is only partially affected by differences in environmental conditions, specifically for the geological, ecological, vegetational and pedological aspects. The soil in the two examined sites is rich in macro- and less rich in microelements, with differences in types and percentages. In particular, the microelements have similar concentrations, except for Zn and Mn. Among the total of 36 compounds detected and identified in the EOs extracted from the analyzed samples at both sites (Difesa di Malta—Sc1; San Basilio—Sc2), the most abundant compounds, often with strong differences between sites, were α-Terpineol, trans-3-caren-2-ol, and Borneol. At Difesa di Malta (Sc1), in more disturbed conditions and in an agriculture context, with a low plant biodiversity referable to an unrepresentative aspect of *Phagnalo saxatilii-Saturejetum cuneifoliae*, to a significant decrease in trans-3-caren-2-ol and borneol, there is a significant increase in α-Terpineol, compared to the better naturalness of San Basilio. Clearly, this means that the environment may play a strong role.

The composition of the *S. cuneifolia* EOs found in the present research does not match with those of previous studies conducted in the Mediterranean area [47,67,69], which showed different chemotypes also among themselves. In Croatia, the chemical composition shows a high percentage of carvacrol (17.7%), γ-terpinene (14.8%), p-cymene (9.8%), linalool (6.6%) and limonene (6.2%) [94]. In another work conducted in Salento (Apulia), but in different areas than the present one, there was a prevalence of linalool (9.6–32.7%), borneol (12.9– 24.0%) and α-pinene (9.5–11.7%) [70]. Recently, in Montenegro, it has been observed that the main constituents were linalool (20.3%), trans-(E)- caryophyllene (6.1%), germacrene D (5.8%), nerolidol (5.2%) and spathulenol (5.0%) [47]. Again, this confirms an environmental effect.

The climatic and geological factors, soil properties, and topographical and vegetational characteristics do not seem to play a crucial role on the chemical composition in EOs of *Thymus spinulosus* in Apulia, as observed in previous studies conducted in a Sicilian population where the environmental conditions and the different soil types affect the composition of EOs. The results of research conducted in Apulia indicate a total of 35 compounds for the three investigated sites, with the highest presence in all samples of thymol, p-Cymene and β-Ocimene, always with comparable percentages, unlike the Sicilian populations, with 62 components in four sites, geographically close and located at comparable altitudes, with variable values of the main components: myrcene (1.0–15.7%), limonene (2.3–13.2%), γ-muurolene (7.3–15.9%), caryophyllene (8.3–11.1%), and germacra-1(10),4-dien-6-ol (1.3–11.3%).

Our results may highlight the following responses in the chemical composition of EOs of *Th. spinulosus*: (a) the Sicilian and Apulian populations are two different chemotypes and the Apulian ones are more stable in composition, regardless of environmental, climatic and vegetational factors; (b) latitude and altitude could be diagnostic factors, even if it would be interesting to complete the picture with other Italian populations such as those of Calabria; (c) it is confirmed that the genus *Th.* has several chemotypes. Furthermore, some studies underscore that other factors are also implicated in determining the EO composition of *Th.*, such as genetic and reproductive characteristics [94], the ecological functions of the EOs,

e.g., protection from herbivores, interaction with microorganisms in the decomposition process, patterns of vegetation through allelopathic action [95], the physiological stage of the plant [87,96], and the negative impacts of herbicides [97].

Finally, adopting specific measures for selective management of the shrub vegetation is recommended, including those identified here, with a high concentration of species rich in EO, especially in the Mediterranean climate, in order to contribute to a reduction in the risk of fire [98] and safeguard the endemic and threatened flora, some of them still without studies of their chemical properties.

As general remarks, for both studied species, we may stress, once more, that the high bioactivity of *Th. spinulosus* indicates it has promising potential in organic agriculture, since it may provide thymol as a possible natural agent against phytopathogenic microorganisms; the richness of α-Pinene (insecticidal) in *S. cuneifolia* also makes the species a potential plant as natural pesticide in organic agriculture.

**Author Contributions:** Conceptualization, E.V.P.; methodology, E.V.P., F.V. and G.N.M.; software, E.V.P., F.V., G.N.M. and A.T.; validation, E.V.P.; formal analysis, E.V.P., F.V., G.N.M. and A.T.; investigation, E.V.P., F.V. and G.N.M.; data curation, E.V.P., F.V., G.N.M. and A.T.; writing—original draft preparation, E.V.P.; writing—review and editing, E.V.P., F.V., G.N.M., S.J. and A.T.; visualization, E.V.P., F.V., G.N.M., S.J. and A.T.; supervision, E.V.P.; pictures, E.V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research did not received funding, but it was performed within the master's program in "Mediterranean Organic Agriculture" assigned to Shaima Jallali at the CIHEAM, Mediterranean Agronomic Institute of Bari.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors wish to thank Lina Al Bitar (CIHEAM) and Ramez Mohamad (CIHEAM) for the administrative support in the management of the Master. The authors are also grateful to the Botanical Garden Museum of the University of Bari, especially Rocco Mariani, the curator of the *Herbarium Horti Botanici Barensis* (BI), Gaetano Pazienza for scanning the specimens of *Satureja cuneifolia* and *Thymus spinulosus*, and Luigi Forte for having authorized their reproduction.

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

#### **References**


## *Article* **Influence of Habitat Types on Diversity and Species Composition of Urban Flora—A Case Study in Serbia**

**Milan Gliši´c 1,2,\*, Ksenija Jakovljevi´c 2, Dmitar Lakuši´c 2, Jasmina Šinžar-Sekuli´c 2, Snežana Vukojiˇci´c 2, Milena Tabaševi´c <sup>2</sup> and Slobodan Jovanovi´c <sup>2</sup>**


**Abstract:** The aim of this study was to investigate the floristic composition and diversity of seven urban habitat types in 24 Serbian cities with different climatic affiliation. In each of the 24 cities, we selected 1 ha plots representing a habitat from one of the following groups: square, boulevard, residential area with compact and with open building pattern, city park, and sites with early and mid-succession vegetation stages. All vascular plant species that occur spontaneously in these plots were observed. Data on the main climatic characteristics were collected for each plot, and data on the life forms were obtained for each species recorded. Diagnostic species were identified for each habitat type analyzed, and alpha, beta and gamma diversity were calculated. A total of 674 taxa were recorded in the studied area. Significant differences were observed in habitats by diagnostic species and by life form representation. The lowest alpha and gamma diversity and the dominance of therophytes were observed in habitat types with intensive anthropogenic impact, whereas the highest number was recorded in mid-successional sites and residential areas with a compact building pattern. The analysis showed that habitat type influences species composition much more than climate.

**Keywords:** species composition; richness; urban areas; anthropogenic impact; climate

#### **1. Introduction**

Human activities are an inseparable part of urban area and play a leading role in modifying its ecological characteristics, forming similar conditions in diverse, often remote areas. Hence, the similar urban habitats are found in the vast majority of cities, even in areas in different biogeographical regions, with different macroclimatic characteristics [1,2]. Large-scale introduction of species with the cosmopolitan type of distribution, sometimes associated with a decline in native species, may lead to a further increasing similarity in species composition between regions [3,4]. Additional homogenization is caused by the presence of invasive species, primarily archeophytes, whereas neophytes mainly lead to opposite effects [5,6]. However, both groups of invasive species have been shown to contribute to an increase in the richness of plant species in urban habitats [5,7]. Namely, according to Pyšek [8], archaeophytes and neophytes account for 15% and 25%, respectively, of the urban flora in Central Europe, although the negative effects of alien species on native diversity have also been observed [9]. Additionally, urban areas are very heterogeneous, so this also contributes to a larger number of species in cities [10,11]. This heterogeneity, caused by different disturbance regimes, induced differences in species composition [12]. Simultaneously with the certain similarities and the number of generalist species found both in and outside the cities, pronounced differences can be observed compared to the surroundings, wherefore the cities can be regarded as a kind of ecological island [6].

**Citation:** Gliši´c, M.; Jakovljevi´c, K.; Lakuši´c, D.; Šinžar-Sekuli´c, J.; Vukojiˇci´c, S.; Tabaševi´c, M.; Jovanovi´c, S. Influence of Habitat Types on Diversity and Species Composition of Urban Flora—A Case Study in Serbia. *Plants* **2021**, *10*, 2572. https:// doi.org/10.3390/plants10122572

Academic Editor: Robert Philipp Wagensommer

Received: 26 October 2021 Accepted: 13 November 2021 Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

There are many characteristics of urban habitats that distinguish them from natural ones: higher levels of disturbance [13,14], herbicide use [15], air and soil pollution [16,17], nutrient enrichment [5,18], higher temperatures due to the urban heat island effect [19,20], higher input of alien species propagules [21], etc.

There are a number of reasons why urban flora research is attracting increasing attention: most of Europe's human population lives in urban areas [22]; the plants of urban habitats contribute to ecosystem services and affect the citizens' well-being [15,23]; urban areas can be centers for the spread of allochthonous plant species to neighboring territories [8,24]; and urban habitats can serve as a refugia for plant species, even those considered rare or endangered [11]. The number of studies dealing with urban flora has increased considerably in recent decades [25]. Comparative studies of urban flora in large areas and the implementation of standardized sampling protocols made a great contribution to the understanding of the distribution and ecology of plants in central European cities [5,6,12,26–29], indicating striking differences between urban habitats in terms of plant species diversity, induced by the types of urban habitats, climate, and specific spatial patterns [26,30,31]. However, the results from Central Europe are hard to generalize to the whole of Europe. For example, diversity of plant species in urban habitats in Southern Europe is greater compared to the parts of Central Europe with a different climate [30]. In addition, the urban flora of Southeastern Europe is poorly studied compared to other parts of Europe, with previous studies often focusing on individual cities or specific urban habitats within them [25]. Comparative and comprehensive studies of urban flora in SE Europe are particularly rare and are not based on standardized sampling methods, rather on a comparison of existing data [32]. For the above reasons, large-scale comparisons and generalizations of features and trends in the urban flora of Southeastern Europe are still lacking.

The aim of this study was therefore to investigate the floristic composition and diversity of urban habitats in Serbia using the standardized protocol established by Lososová et al. [26], to obtain comparable results. Cities in Serbia are a good model for such a study, as climatic differences between cities are considerable. Hence, our aim was also to investigate (1) how urban habitat types affect floristic composition and diversity and (2) which factor has a stronger effect on floristic composition, habitat types, or climatic characteristics.

#### **2. Results**

Significant differences in species composition were observed among selected habitat types. Diagnostic species for each habitat type are listed in Table 1. The highest average value of Φ (0.13) and the highest number of diagnostic species (67) were calculated for mid-successional sites (m). This habitat type was dominated by deciduous shrubs and trees (e.g., *Prunus spinosa*, *Cornus sanguinea*, and *Juglans regia*), but also by a large number of herbaceous perennial species (e.g., *Hypericum perforatum*, *Dipsacus laciniatus*, *Rumex patientia*, *Agrimonia eupatoria*) and grasses (e.g., *Calamagrostis epigejos*, *Poa trivialis*, *Holcus lanatus*, *Bromus hordeaceus*). Early successional sites and residential areas with compact building patterns were also characterized by a relatively high average value of Φ (0.10, both) and a relatively high number of diagnostic species (39, both). The diagnostic species with the highest Φ value (>0.50) in early successional sites were *Papaver rhoeas*, *Bromus hordeaceus*, *Polygonum lapathifolium*, and *Vicia cracca*. Residential areas with compact building patterns were mostly characterized by ornamental species, which spread from neighboring gardens (e.g., *Campsis radicans*, *Kerria japonica*, *Antirrhinum majus*, *Ipomoea purpurea*, *Rudbeckia hirta*). In residential areas with an open building pattern, 22 diagnostic species were observed, with domination of ornamental species (e.g., *Mirabilis jalapa*, *Commelina communis*, *Hibiscus syriacus*) and an average Φ value of 0.7. Only five diagnostic species were determined for city parks (p) and the species with the highest Φ value was *Quercus robur*, which often occurs as seedlings. Boulevards (b) hosted a large group of different plants, but only four species were determined as a diagnostic. By comparison, no positive correlation with city

squares (s) was found for any of the plant species, and the lowest average value of Φ was obtained (0.01) for this habitat type.

**Table 1.** Diagnostic plant species for studied urban habitat types in Serbia.


Taxa are listed by decreasing values of Φ (in parenthesis). Only taxa with Φ > 0.3 are shown. The taxa with Φ > 0.5 are marked bold.

The habitat types also differed significantly in the representation of plant life forms (Figure 1). The largest participation of hemicryptophytes, as most represented compared to the other life forms, was found in early and mid-successional sites, and the lowest at squares and boulevards. By comparison, therophytes had the largest share in squares and

the lowest in mid-successional sites. Parks and residential areas, both compact and open, were characterized by a somewhat higher proportion of phanerophytes and chamaephytes compared to other habitat types. Geophytes were generally the least represented, especially in mid-successional sites.

**Figure 1.** Proportion of hemicryptophytes, therophytes, phanerophytes/chamaephytes, and geophytes in particular urban habitat types. X-axis abbreviations: b—boulevard, e—early successional sites, m—mid-successional sites, p—city park, c—residential area (compact building pattern), o—residential area (open building pattern), s—historical city square. Homogeneous groups of urban habitat types are denoted by the same letters (*p* < 0.01).

The ordination diagram obtained by PCA indicated a predominant grouping of plots belonging to the same habitat type, with a certain overlapping between some of them (Figure 2). The city square (s) distinguished most clearly from the others, with no overlap of this habitat type with any other habitat type analyzed. The most distant from the other habitat types were the plots representing mid-successional and early successional sites. These two habitat types formed a well-separated group, with minor overlapping between each other. The differences were much less noticeable between other habitat types, particularly between city parks (p) and boulevards (b) and between the two types of residential areas (c and o).

The results of the RDA analysis indicated that species composition was affected by both habitat type and climatic variables, with the greater influence of habitat type (9.7%) in relation to climatic characteristics (2.3%), whereas shared variance was not observed (Table 3).

The total number of taxa recorded on 164 plots in the 24 investigated cities on the territory of Serbia was 674, with an average of 105 taxa observed per plot. The highest number of taxa was recorded in the residential areas with compact building pattern (c) in Niš (147), Kruševac (146), and Loznica (142), whereas the lowest number of registered taxa was observed in the city squares (s) of Smederevo and Panˇcevo (38) and Sremska Mitrovica (42).

There were considerable differences in alpha, beta, and gamma diversity between the habitat types studied. The lowest alpha diversity was found in the city squares (s; 53 plant taxa) and boulevards (b; 88 taxa), whereas the highest values were recorded in residential areas with compact building patterns (c; 128 plant taxa) and at the mid-successional sites (m; 127 plant taxa). Similarly, the lowest gamma diversity was found in the city squares (s; 244 plant taxa) and boulevards (b; 297 taxa), whereas the highest values were found at the mid-successional sites (m; 435 plant taxa) and residential areas with compact building pattern (e; 398; Figure 3A). The highest beta diversity was recorded within city squares (s) and the lowest in residential areas with open building patterns (c; Figure 3B).

**Figure 2.** PCA ordination of plots according to plant species composition. Eigenvalues: axis 1, 0.2492; axis 2, 0.0925. Abbreviations: the first two uppercase letters represent the city code (see Table 2), the third lowercase letter represents the habitat type (see Material and Methods). Legend: red triangles boulevards, orange squares—early successional sites, light green squares—mid-successional sites, dark green boxes—parks, purple circles—residential area (compact building pattern), light blue circles—residential area (open building pattern), blue diamonds—squares.



**Table 2.** *Cont*.

N—latitude, E—longitude, Alt—elevation, T—mean annual temperature, ΔT—difference between mean temperature in July and January, P—precipitation.

**Table 3.** Total explained variation, the influence of habitat type, climate, and their shared effect on the plant species composition in the analyzed cities.

**Figure 3.** (**A**) Alpha (box plots) and gamma (numbers above) diversity in particular urban habitat types. (**B**) Beta diversity of plant taxa within urban habitat types. X-axis abbreviations: b—boulevard, e—early successional sites, m—mid-successional sites, p—city park, c—residential area (compact building pattern), o—residential area (open building pattern), s—historical city square. Homogeneous groups of urban habitat types are denoted by the same letters (*p* < 0.01).

A strong positive statistically significant correlation was found between alpha and gamma diversity (r = 0.93, *p* = 0.002). Alpha and beta diversity were strongly negatively and statistically significantly correlated (r = −0.82, *p* = 0.02). Gamma and beta diversity were also negatively correlated, but this correlation was statistically insignificant (r = −0.62, *p* > 0.1).

#### **3. Discussion**

Significant differences in species composition were found between individual types of urban habitats. As observed from the results of ordination analysis and the plots grouping, species composition is primarily influenced by urban habitat type and intensity of human influence, much more than by climatic features. This is consistent with the findings of Lososová et al. [26], who analyzed the urban flora of Central European cities along a gradient of distinctly different biogeographical regions with contrasting climatic characteristics, and with the results of Rebele [1] and Savard et al. [2], which confirm the hypothesis of uniformity of the urban environment.

Based on the ordination analysis of species composition, plots were grouped primarily based on habitat type affiliation. Four groups of plots can be observed: plots representing squares (1), plots representing boulevards and parks (2), plots representing residential areas (3), and plots representing successional sites on the urban peripheries (4). This suggests that certain habitat types have greater similarities in species composition than others. These similarities are the result of a similar character of anthropogenic influence and the location of the plots in the city (center or periphery; Figure 2).

The results of alpha and gamma diversity analysis in the urban habitats of the studied area showed that the lowest number of species was found in city square, which is consistent with the results from Central European cities [6,27]. Sealed and paved surfaces, which dominate in city squares are hostile habitats, and a relatively small number of plant species are adapted to such extreme conditions (trampling, high insolation, drought, etc.). Hence, it is not surprising that therophytes, as a disturbance tolerant life form, have the largest share in city squares (Figure 1). However, although well-adapted to highly influenced urban habitats, therophytes are thought to be more prone to extinction compared to other life forms [33,34]. Intensive human influence and vegetation limited to small patches strongly support the presence of cosmopolitan and even alien species, which negatively affect the native flora and consequently leads to their local extinction [34,35]. The richness of native flora is additionally affected by environmental filters that have to be overcome in order for plants to arrive from the surrounding natural habitats [12]. Furthermore, these habitat types have a significant share of alien and ornamental plants that have been spread by human activities. Due to all this, squares are unique in their species composition compared to the other urban habitat types, as indicated by their clear separation from them in the ordination diagram (Figure 2).

As opposed to a small number of species in city squares in Serbia and in Central European cities [27], the high plant species richness has been observed in the old centers of the Italian Mediterranean cities [30], due to a number of plants growing undisturbed on the ancient walls. However, such microhabitats beyond intense anthropogenic influences are rare in the city squares in Serbia, and chasmophytic flora were represented by only a few species (e.g., *Cymbalaria muralis* and *Sedum* spp.). Simultaneously with the lowest alpha and gamma diversity, urban squares in Serbia had the highest beta diversity, indicating that the differences between specific habitats of this type were greater than in other habitat types. The same was found in Central European cities as a result of a very low number of species [26].

A somewhat higher floristic richness was observed in boulevards. Bearing in mind that boulevards are usually spatially connected to city squares, the human influence is similar to that for squares but of lower intensity, which resulted in somewhat higher species richness. A major contributor to species richness in boulevards is tree lines, especially those with unsealed surfaces around the trees, considering that these microhabitats harbored a wide range of plant species. Due to the presence of planted trees, a significant number of their seedlings were observed in boulevards. This is one of the reasons for the marked similarities in the species composition of boulevards and parks, whose plots are grouped in the ordination diagram (Figure 2). The species detected in this habitat type are generally common in the cities and almost half of them were detected in all the habitat types analyzed, whereas only nine were recorded only in the boulevards, confirming that the species found in the boulevards are those that generally survive very successfully in urban conditions [12]. Similar to the city squares, the boulevard's flora were characterized by a significant participation of therophytes, particularly compared to the other habitat types. However, in addition to therophytes, boulevards were also characterized by perennials adapted to growing in small cracks in concrete (e.g., *Sagina procumbens*) and seedlings of tree species, primarily those deliberately planted (e.g., *Acer pseudoplatanus*). In contrast, in city squares and boulevards, the least represented life form was geophytes (Figure 1). Although bulbs and rhizomes enable them to survive in the hostile environment [36], their numbers were found to be negatively correlated with the number of inhabitants and traffic density [32].

Parks are very specific urban habitats. Although they resemble natural forest habitats (e.g., similar mild-mesoclimate, similar light regime, etc.), with a number of typical forest plants that can be found in them (e.g., *Brachypodium sylvaticum*, *Clematis vitalba*, seedlings of *Quercus robur*, etc.), parks are artificial habitats created by deliberate planting of trees,

often of non-native origin and heavily impacted by human activity (e.g., trampling, regular mowing, planting horticultural annuals and perennials, etc.). In addition to the planted tree species, city parks in Serbia are also characterized by significant participation of spontaneously growing ornamental alien plants, especially typical ruderal ones (e.g., *Amaranthus deflexus* and *Lactuca serriola*, etc.), but also shrubs and tree seedlings (e.g., *Philadelphus coronarius*, *Symphoricarpos albus*, and *Broussonetia papyrifera*). Hence, in addition to phanerophytes as the predominant life form, a significant number of hemicryptophytic species occur in this habitat type (Figure 1). Contrary to almost half of the species detected in the urban parks that were found in all seven groups of the analyzed urban habitat types, a total of 19 species, mainly tree and shrub species, were detected exclusively in this habitat type. However, alpha and gamma diversity in city parks in Serbia was not high, especially compared to residential areas and early and mid-successional sites, primarily because their homogeneity and the lack of specific microhabitats that could harbor different species types. Regardless of this, city parks are recognized as very important habitats for conservation of local biodiversity [37–39].

Residential areas of cities in Serbia, especially those with compact building patterns, harbor a wide range of species, with the highest alpha and significant gamma diversity observed in this type of urban habitat. This high diversity is a consequence of the marked heterogeneity of habitats, considering that residential areas with densely distributed individual housing units represent a mosaic of diverse habitat types: sidewalks and other paved and sealed surfaces, lawns with different mowing regimes, ornamental gardens, urban gardens with various cultivated plants and accompanying weed flora, tree lines, etc., which contribute significantly to the high species richness, especially when compared to squares, boulevards, and urban parks [40]. However, it should be considered that the shape of the sampling plot may also have an influence on increasing species richness and that elongated plots, such as those used to study of this type of habitat, tend to have more species than compact plots of the same size [41–43]. Although our analyses revealed that the highest alpha diversity was found in residential areas with compact building patterns, Godefroid and Koedam [44] indicated a contrasting pattern, i.e., that half open and open areas in Brussels promote species richness, whereas areas with compact structures lead to its reduction. Due to the complexity of these habitat types, residential areas in analyzed cities in Serbia were characterized by different groups of plants, including species commonly found in squares or boulevards (e.g., *Arenaria serpyllifolia* and *Eragrostis minor*), a wide range of grasses (e.g., *Poa* spp. and *Lolium perenne*), escaped ornamental plants (e.g., *Kerria japonica* and *Antirrhinum majus*), juveniles of crop plants (e.g., *Zea mays* and *Solanum lycopersicum*), crop weeds (e.g., *Cynodon dactylon* and *Elymus repens*), and spontaneously growing native and alien trees and shrubs (e.g., *Campsis radicans*, *Syringa vulgaris*, and *Rhus typhina*).

Urban habitat types that stood out in terms of species composition and plant diversity were successional sites, particularly mid-successional ones, abandoned long enough to form species-rich grasslands. Similar findings have been made in several other urban studies [26,44]. Considering that these habitats usually develop on the urban periphery, they are characterized by the absence of strong human impact, with the main difference between them being the duration of the disturbance-free period. Due to recent disturbance, early successional sites are characterized by a vegetation cover dominated by annual plants with ruderal life strategies, capable of rapid colonization of bare ground (e.g., *Bromus tectorum* and *Petrorhagia prolifera*). Additionally, the results of previous studies [35,45] show that compared to the older succession stages, the younger ones are more susceptible to alien species invasion. Mid-successional sites hosted a higher number of plant species than early successional ones and are characterized by the highest alpha diversity compared to the urban habitat types studied. Because of a longer period of non-disturbance, these sites are suitable for inhabiting species with different life strategies. Hence, a lower proportion of therophytes was observed at the mid-successional sites compared to the early successional ones, whereas both habitat types were characterized by a greater participation

of hemicryptophytes and therophytes in relation to the other types analyzed (Figure 1). Additionally, proximity to natural vegetation and openness to the urban surroundings facilitated the influx of native species into these habitats, and various shrub and tree species are common in these habitat types (e.g., *Prunus* spp., *Ulmus* spp., and *Juglans regia*). In general, mid-successional sites are characterized by a higher share of native flora than other urban habitat types, particularly city squares and boulevards [5,35]. Despite the very different physiognomy, early and mid-successional sites have a very similar species composition. Many of the same species of shrubs and trees are found in both types of successional sites, but in the early successional sites they appear as seedlings and much smaller individuals. Additionally, their similarity is also contributed by the fact that due to the complex human influence of varying intensity, some patches of older vegetation were found in the early successional sites analyzed because it was hard to find a completely uniform 1 ha area.

Although the results of this study indicated that local site conditions are the predominant factors determining plant species composition, the results of similar studies of urban flora in Italy contradict this hypothesis because climate was found to be the main factor determining species composition, most likely due to the strong climatic gradient from the north to the south of the Apennine Peninsula [30]. Despite the fact that the climate in Serbia is very diverse, ranging from continental to sub-Mediterranean and mountainous, the differences are not pronounced enough to point to climate as the main factor determining the diversification of flora in urban habitats, which is also true for the Central European urban flora [26]. However, the effect of the climate on the composition of plant species in urban habitats in Serbia should not be neglected. This applies particularly to Sjenica, which is located at an altitude of 1026 m and is characterized by a humid mountain alpine climate (Table 2), as the ordination analysis indicated a grouping of several plots of this city. The climatic differences between the other cities analyzed are less noticeable, resulting in a greater similarity in the ordination diagram (Figure 2). Apart from the climate, the composition of species in urban habitats may be influenced by other factors, such as the geographical location, and the structure of the city and its historical features. However, as their effects can be significant, especially in large-scale studies [27], and almost negligible in studies focusing on cities from a smaller geographical area, these factors are usually neglected in the analyses.

#### **4. Materials and Methods**

#### *4.1. Data Sampling*

Investigation of urban flora was carried out in 24 cities in Serbia (Table 2). Cities were selected to represent all major climate types and subtypes in Serbia. An additional condition that cities had to meet was the existence of the preselected typical urban habitat types. According to Stevanovi´c and Stevanovi´c [46], the following climate types and subtypes are represented in the territory of Serbia (climate type and subtype designations, according to Walter and Leith [47] and Horvat et al. [48], are given in parenthesis): transitional submediterranean-Aegean subcontinental (IV6), semi-arid continental Pannonian (VII), semi-humid continental Danubian (VII), transitional subcontinental-semiarid continental (VI3b/VII), semi-arid temperate continental (subcontinental)—central-southeastern Balkan or Moesian (IV3), humid temperate-continental—west Balkan or Illyrian (IV2b), and humid mountain alpine (XI).

To make comparable samples, floristic data were collected using a standardized protocol, which has already been used in similar studies for Central European cities [5,6,12,26–29]. In each of 24 selected cities in Serbia, we recorded the plant composition at seven specific sites of 1 ha in size. Each of the seven sites represents a different type of urban habitat:


Adequate sites (habitats) were selected using maps and satellite images in Google Earth. The 1 ha plots were selected within each habitat type and all vascular plants were recorded within them. This included seedlings from spontaneously grown planted trees and garden plants. However, intentionally planted individuals were omitted. In residential areas with a compact building pattern (c), a different approach was applied due to limited access to private gardens. In these cases, instead of an area of 1 ha, transects along a street 500 m long were analyzed, recording all species found in accessible areas, in addition to all those that could be seen from the street inside private yards. The research was conducted in 2015–2019 in the period from June to the end of August, to avoid spring and autumn, i.e., plants with significant variations in phenology. The nomenclature of species corresponding to the diagnostic species of classes of plant communities dominated by vascular plants follows Electronic Appendix S6 (EVC1) of Vegetation of Europe [49], and for other species follows the nomenclature of Flora Europaea [50].

#### *4.2. Data Analysis*

The composition of plant species recorded within individual habitat types was shown in synoptic tables. To determine diagnostic species for particular habitat types, the phi coefficient of association (Φ) was used as a statistical measure of the concentration of occurrence of species in particular habitat types [51]. Diagnostic species for a particular habitat type were defined as species that preferentially occur in that habitat type. Fisher's exact test (*p* < 0.05) was used to assess the statistical significance of the species-habitat association, quantified by Φ, as shown in the following equation [52]:

$$\Phi = \frac{\mathbf{N} \times \mathbf{n}\_{\mathrm{P}} - \mathbf{n} \times \mathbf{N}\_{\mathrm{P}}}{\sqrt{\mathbf{n} \times \mathbf{N}\_{\mathrm{P}} \times (\mathbf{N} - \mathbf{n}) \times (\mathbf{N} - \mathbf{N}\_{\mathrm{P}})}}$$

where N is the number of all sites in the data set, Np is the number of sites in the particular habitat type, n is the number of occurrences of the species in the data set, and np is the number of occurrences of the species in the particular habitat type. Diagnostic species were considered to be those that had a statistically significant species–habitat association and Φ > 0.30. Synoptic tables and calculations of phi coefficient were carried out in the JUICE program [53].

All recorded species were categorized into four different groups according to their life form: geophytes, phanerophytes and chamaephytes (trees and shrubs), therophytes, and hemicryptophytes [54]. To compare the differences in the frequency of individual life forms in relation to urban habitat types, ANOVA [55] was used.

According to detrended correspondence analysis (DCA), the length of the gradient in species composition was 2.27 SD units. Therefore, to assess the general variation patterns in composition of plant species among urban habitat types, unconstrained and constrained linear ordination methods were used (PCA and RDA, respectively). Analysis and visualization of the PCA and RDA diagrams were performed using CANOCO 5.12 program [56].

To distinguish the influences of climatic variables and urban habitat types on plant composition, redundancy variation partitioning for RDA was carried out [57]. Two groups of variables were collected for each site: (1) habitat type, which was given as a categorical variable with seven expressions, and (2) three climatic variables: mean annual temperature, annual temperature range, i.e., the difference between mean temperatures in July and January, and annual precipitation total. The climate variables were taken from the WorldClim dataset. The significance of the influence of climatic variables and urban habitat types was tested by Monte Carlo permutations (999 permutations). These calculations were performed in the CANOCO 5.12 program [56].

Alpha and gamma diversity were used to indicate differences in plant species richness between habitat types. Alpha diversity was defined as the average number of taxa recorded per plot in each habitat type, and gamma diversity was estimated as the total number of taxa observed in all plots belonging to a particular habitat type. To determine beta diversity, we calculated an index of beta diversity: S/a-1, where S is the total number of taxa, and a is the average number of taxa per plot. The calculations of alpha, gamma, and beta diversity and the visualization of their differences between habitat types were carried out in the "vegan" package of the R programming language [58].

#### **5. Conclusions**

The results of the study of urban flora of Serbia indicate significant floristic richness of the investigated areas. Pronounced differences were found between the analyzed habitat types, both in terms of diversity, species composition, and dominant life forms. It was shown that the species composition and dominance of particular life forms are directly related to the intensity of anthropogenic influence. The lowest alpha and gamma diversity was found in city squares and boulevards, with a dominance of therophytes, a life form related to disturbed habitats, whereas habitats in the urban periphery under less pronounced human influence and heterogeneous residential urban areas, are characterized by significant floristic richness and much more uniform distribution of life forms. Comparison of the influence of different factors on the composition of urban flora shows that urban flora in the analyzed cities is more influenced by the type of habitat than by climatic features. Our results indicated that the most important ways to increase plant diversity in cities are the following: (1) allowing natural succession and reducing the intensity of the anthropogenic factor in certain parts of the city; (2) providing greater heterogeneity of urban areas, and forming specific urban microhabitats that allow the survival of those species that are not typical urbanophiles.

**Author Contributions:** Conceptualization, M.G., K.J., J.Š.-S., D.L. and S.J.; writing—original draft preparation, M.G., K.J., J.Š.-S. and S.J.; methodology, M.G., J.Š.-S., D.L. and K.J.; software, M.G. and J.Š.-S.; investigation, M.G. and M.T.; validation, M.G., D.L. and S.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Ministry of Education, Science and Technological Development of the Republic of Serbia, grant number 451-03-9/2021-14/ 200178.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are not publicly available due to the database is part of the Ph.D. thesis. The complete database will be available only after the publication of the dissertation in 2022. Also, we will use the same database for other articles we plan to publish. For these reasons, we believe that the complete database should not be publicly available at this moment.

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

#### **References**


### *Article* **Orchids of Azerbaijani Cemeteries**

**Attila Molnár V. 1,\*, Viktor Löki 2, Marc Verbeeck <sup>3</sup> and Kristóf Süveges <sup>1</sup>**


**Abstract:** In order to explore their orchid flora, we performed surveys of 96 Azerbaijani burial places in 2018 and 2019. Altogether, 28 orchid taxa were found in 37 visited cemeteries. In the orchid diversity a remarkable pattern was observed: geographic latitude was significantly and positively related to the number of taxa and number of individuals. The most widespread and abundant orchids in Azerbaijani graveyards were *Anacamptis pyramidalis* and *A. papilionacea* (found in 23 and 8 cemeteries, respectively). Azerbaijani cemeteries can be important refuges for rare and threatened orchids, e.g., *Himantoglossum formosum* (three cemeteries), *Ophrys sphegodes* subsp. *mammosa* (eight), *Orchis adenocheila* (two), *O. punctulata* (three), *O. stevenii* (one) and *Steveniella satyrioides* (one). *Epipactis turcica*, detected in a single locality, was previously unknown to the flora of Azerbaijan. Additionally, we documented orchid tuber (salep) collection in two cemeteries.

**Keywords:** anthropogenic habitats; Caucasus; *Himantoglossum formosum*; human-made habitats; Orchidaceae; salep harvesting; Transcaucasia

#### **1. Introduction**

The Earth's surface has changed dramatically in recent centuries, with human activities serve as a leading cause of the drastic reduction in the area of natural habitats [1,2]. In parallel with the degradation and fragmentation of natural environments throughout the world, isolated natural habitat patches as remnants of the original wildlife have been revalued [3]. Anthropogenically influenced habitats now occupy a significant part of the Earth's surface and expand rapidly [4]. In order to conserve the remaining biodiversity, it is of the utmost importance to identify and protect the remaining habitats with a high conservation value, to develop a sustainable habitat management practice, and to plan future developments in the light of nature conservation priorities [5].

Recently, conservation professionals have recognized that some of the anthropogenically influenced or even human-made habitats, such as abandoned mines and industrial sites [6–8], road verges [9–11], tree plantations [12–14], river dikes [15], burial mounds [16], and urban habitats [17,18], play significant roles in conserving biodiversity. During the last decades, it has become increasingly evident that cemeteries also play an important role in maintaining biodiversity [19]. Although the orchid flora of cemeteries is globally rather poorly known, occurrences of orchids were published from Australian, Asian, and European burial places [20]. Based on previous knowledge on the occurrence and diversity of orchids in Turkish [21–24], Albanian [25] and central European [26] burial grounds, we predicted potential conservational importance of traditional Caucasian cemeteries. One of the main goals of our study was to search for *Himantoglossum formosum*, the rarest and perhaps the least known orchid of the Caucasian region [27]. During the 180 years after its description [28], almost nothing was known about the species [29], and it was re-discovered in 1994 [30]. According to recent studies, this is an ancient, phylogenetically isolated [31]

**Citation:** Molnár V., A.; Löki, V.; Verbeeck, M.; Süveges, K. Orchids of Azerbaijani Cemeteries. *Plants* **2021**, *10*, 2779. https://doi.org/10.3390/ plants10122779

Academic Editor: Robert Philipp Wagensommer

Received: 30 November 2021 Accepted: 11 December 2021 Published: 16 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and morphologically well separated [32] bona fide species. It is listed as Vulnerable (Rare) in the IUCN Red List of Threatened Plants [33].

The aims of this paper were to survey Azerbaijani cemeteries as orchid habitats, and to test which geographic factors influence the prevalence of orchids in the surveyed cemeteries.

#### **2. Materials and Methods**

We studied burial grounds (Azerbaijani: m@zarlıq, hereafter cemeteries) regardless of their spatial dimension, position within settlements, or presence of built facilities. We surveyed 96 Azerbaijani cemeteries (Figure 1, Table A1) during 2018 (17–30 May by Molnár V., Löki, Mizsei and Süveges, and 28 June–4 July by Molnár V. and Szabó) and 2019 (29 April–6 May by Verbeeck, Duijnhouwer, Segers and Bobocea) and (31 May–6 June by Verbeeck, Duijnhouwer and Bradeanu). Most cemeteries were visited only once (90 and 3 cemeteries in May 2018 and in April 2019, respectively), but three cemeteries were visited in both years. All orchid taxa and the number of individuals were counted or estimated in the whole area of each visited cemetery. Species were identified based on the comprehensive book of Kuehn et al. [34]. Authors of plant names were listed in Table 1. The geocoordinates and the elevation of the visited cemeteries were determined using a Garmin eTrex Legend handheld GPS device and recorded in WGS84 format. During field trips, particular attention was devoted to documenting salep collection activity in cemeteries.

**Figure 1.** Number of orchid taxa in the cemeteries surveyed.


**Table 1.** Orchid taxa recorded in Azerbaijani cemeteries.

To understand the role of geographic factors in determining variation in taxon richness and abundance of orchids across Azerbaijan, we built statistical models with either of these variables as dependent variables, and latitude, longitude and altitude as explanatory variables. Both the number of individuals and the number of taxa had Poisson distributions, but due to the overdispersion in these variables, we used generalized linear model (GLMs) with quasi-Poisson distribution. All models were built in the R statistical environment [35].

#### **3. Results**

Numbering (ID), geographic location, and altitude above see level of the cemeteries visited, together with lists of the orchid taxa found in each one, are given in Table A1. In total, 28 orchid taxa were found, and considerable differences can be observed in the number of individuals and frequency of each taxon (Table 1), as well as in orchid species richness and abundance of each cemetery (Table 2).

**Table 2.** Descriptive statistics orchid flora of Azerbaijani cemeteries.


Each taxon was found total in 1–24 cemeteries (mean ± SD = 3.2 ± 4.5), with the number of individuals varying from 1 to 1902 (mean ± SD = 150 ± 374). The most widespread and abundant species was *Anacamptis pyramidalis* (Figure 2A). The number of taxa detected in only one graveyard was 15, whereas four species were found in more than five cemeteries. The highest number of taxa in a given cemetery was 9. In most cases only one taxon (18 cemeteries (15%)) or two taxa (11 cemeteries (9.4%)) occurred. Cemeteries that serve as habitats for five or more taxa were extremely rare (4 (3.4%)). The most orchid-rich cemeteries were found near Lerik (AZ-16, 9 species) A ˘gab@yli (AZ-52, 8 species), Nohurqishlaq (AZ-93, Figure 2B, 8 species), and DashliJalgan (AZ-90, 5 species).

**Figure 2.** Orchids in Azerbaijani cemeteries. (**A**) *Anacamptis pyramidalis* population in the cemetery of @ngixaran (AZ–61). (**B**) Cemetery of Nohurqishlaq (AZ-93), habitat of *Orchis punctulata*, *O. stevenii* and their hybrids (*Orchis* ×*chabalensis*). (**C**) Viable population of *Himantoglossum formosum* was found on a few tens of square meter of refuge under some old oak trees in cemetery of Zizik (AZ-74). (**D**) Spurs of salep harvesting in the cemetery of A ˘gab@yli (AZ-52). (**E**) Effect of fencing around cemetery against grazing: plant cover is considerable lower outside (left) than inside (right, with flowering individuals of *Anacamptis pyramidalis*) of cemetery of Zurnabad (AZ-32). (**F**) Inflorescence of *Steveniella satyrioides*. (**G**) Occurrence of *Epipactis turcica* was formerly unknown from Azerbaijan (Tengealti, AZ-85). (**H**) A very localized and rare endemic species, *Himantoglossum formosum* in cemetery of Zizik (AZ-74). Photo credit: **A**, **C**, **D**, **G** and **H** by A. Molnár V.; **B** and **F** by M. Verbeeck; **E** by V. Löki.

The harvest of orchid tubers ("salep") was observed in two cemeteries during 2018. In A ˘gab@yli cemetery (AZ-52, Figure 2D) three species (*Anacamptis papilionacea, Orchis adenocheila, O. simia*), and in Dashli Jalgan cemetery (AZ-90) five species (*Anacamptis collina,* *A. papilionacea, Ophrys sphegodes* subsp. *mammosa, Orchis simia, Neotinea tridentata*), were collected. Both of these localities host notable orchid populations with eight and five species, respectively.

The number of orchid taxa and individuals found in Azerbaijani cemeteries was significantly positively related to latitude (Tables 3 and 4, respectively), but not to longitude and altitude. When non-significant predictors were removed from the model in a stepwise manner (based on the largest *p*-values), only latitude remained in the final model as a significant predictor of orchid species richness and abundance.

**Table 3.** Effect of geographic location on number of orchid taxa per cemetery. Parameter estimates, their standard errors (SE), associated t-values (t) and significance levels (p) are presented.


**Table 4.** Effect of geographic location on number of orchid individuals in Azerbaijani cemeteries.


#### **4. Discussion**

During our work, it has been proved that Muslim Azerbaijani cemeteries host significant orchid populations. The key conservation importance of Azerbaijani cemeteries can be explained by two facts: (1) Religious privileges protected these sacred sites and their natural values, because they have largely been exempt from forest and agricultural utilization ever since; and (2) the mostly fenced area of cemeteries provide protection against excessive grazing (Figure 2E).

Azerbaijani cemeteries provide shelters for several valuable populations of rare and threatened orchids. From a conservation point of view, one of the most valuable species is the Eastern Caucasian endemic *Himantoglossum formosum* (Figure 2H), which was found in three of the visited cemeteries (Zizik, AZ-74, Figure 2C; Yasab, AZ-78; Piral, AZ-79). Viable populations of the rare *Orchis adenocheila* were found in two cemeteries (Lerik, AZ-16; A ˘gab@yli, AZ-52). The occurrence of *Steveniella satyrioides* was detected in cemetery of Lerik (Lerik, AZ-16, Figure 2F). The occurrence of *Epipactis turcica* (Figure 2G) was also found near Tengealti (AZ-85); this taxon was formerly unreported in Azerbaijan.

The long-term survival of these orchid populations in cemeteries strongly depends on long-established, sustainable management practices and traditional burial habits [22,36]. Establishment of graves (especially modern graves covered by marble or concrete tombstones) on the most valuable parts of these cemeteries is expressly undesirable from a conservation perspective, as well as the use of herbicides or electric trimmers. However, mowing or moderate grazing of grassy areas around the burial ground is preferred and encouraged for a more efficient conservation of the local biodiversity and valuable flora elements. Based on their diverse and abundant orchid community in some of the visited cemeteries, we strongly recommend the local councils and the nature protection authorities to protect certain burial places, especially near Lerik (AZ-16), A ˘gab@yli (AZ-52), DashliJalgan (AZ-90), Nohurqishlaq (AZ-93), and Nugadi (AZ-92).

A special threatening factor of tuberous orchids, namely the harvest of their tubers (making salep for culinary purposes [37]) was observed in Azerbaijani cemeteries. On

the one hand, the right of local human communities to continue using traditional natural resources is unquestionable and seems also sustainable [38,39]. On the other hand, the effects of tuber collection on populations of frequent and widespread orchids is little known, while the sustainability of salep harvesting is at least controversial [40–47]. However, destroying the rarest taxa (*Himantoglossum formosum, Orchis adenocheila*) should definitely be avoided.

**Author Contributions:** Conceptualization: A.M.V.; formal analysis: K.S.; investigation: A.M.V., V.L., M.V. and K.S.; writing—original draft preparation: A.M.V.; writing—review and editing: V.L., M.V. and K.S.; visualization: A.M.V., V.L. and M.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Research, Development and Innovation Office of Hungary (grant number NKFI-OTKA K132573).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data analyzed in this study are available in Appendix A.

**Acknowledgments:** The authors are grateful to Edvárd Mizsei, Éva Szabó, Luc Segers, Roel Duijnhouwer, Alin Bradeanu and Mihai Bobocea for their assistance during the field work, to Attila Takács for editing the map of surveyed cemeteries and to the anonymous reviewers for their valuable suggestions and useful recommendations. We would like to express our gratitude to C. A. J. Kreutz (The Netherlands) for identification of *Epipactis turcica*, and to Orsolya Vincze for her linguistic corrections of the earlier version of our manuscript.

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

#### **Appendix A**

**Table A1.** Numbering (ID), geographic location, altitude, year of observation and orchid taxa of the 96 cemeteries studied in Azerbaijan. A dash "–" indicates that no orchid taxa were recorded. Generic name abbreviations: *A.—Anacamptis, C.— Cephalanthera, D.—Dactylorhiza, E.—Epipactis, H.—Himantoglossum, L.—Limodorum, O.—Orchis, Op.—Ophrys, S.—Steveniella*.



#### **Table A1.** *Cont.*


#### **Table A1.** *Cont.*

#### **References**

