Monte Argentario is a coastal promontory located in southwestern Tuscany, central Italy (42.399339°N, 11.154622°E) (Fig. 1). Although orographically isolated from the mainland, it is connected to the Tuscan coast by two sandy tombolos and an internal lagoon system (Lazzarotto et al. 1964). The promontory covers approximately 60 km² and is characterized by steep rocky cliffs and rugged topography, reaching 635 m a.s.l. at Punta Telegrafo.

Geologically, the area exhibits a complex assemblage of limestone and dolostone, metasandstone, acid metavolcanic rock, calcschist, and metamorphic ophiolite (Carmignani and Lazzarotto 2004), while the climate is Mediterranean (Csa type: hot-summer Mediterranean climate; Köppen 1936), with mild winters and hot, dry summers (Arrigoni and Di Tommaso 1997).

Bioclimatically, the area falls within the Mediterranean macrobioclimate (pluviseasonal-oceanic variant), spanning a gradient from the lower mesomediterranean thermotype at coastal elevations to the upper mesomediterranean thermotype towards the summit (Pesaresi et al. 2017), within the Mediterranean Biogeographic Region (European Environment Agency 2017).

The vascular flora of Monte Argentario includes approximately 1200 taxa. The biological spectrum (Raunkiaer 1934) shows a predominance of therophytes (40%), followed by hemicryptophytes (30%), and geophytes (15%). This pattern, together with the prevalence of Mediterranean, East-Mediterranean, and Atlantic-Mediterranean chorotypes confirms the strictly Mediterranean character of the area (Baldini 1995).

The structure of the vegetation is a dynamic mosaic of sclerophyllous evergreen forests and their degradation stages, including Mediterranean shrubland, garrigues, and discontinuous annual grasslands. This mosaic reflects long-term anthropogenic influence (agriculture, grazing, recurrent fires, and land abandonment) superimposed on natural successional processes.

Phytosociologically, these communities can be framed within the classes Quercetea ilicis Br.-Bl. ex A. Bolòs y Vayreda 1950, Querco-Fagetea Br.-Bl. and Vlieger in Vlieger 1937, and Festuco-Brometea Br.-Bl. and Tx. ex Br.-Bl. 1949 (Arrigoni and Di Tommaso 1997).

Thirteen habitat types listed in Annex I of the EU Habitats Directive (92/43/EEC) have been identified in the area (Viciani et al. 2018).

The most representative habitats are:

Habitats 5330 and 6220* frequently form intricate mosaics representing dynamic successional or degradation stages. Abandoned olive groves (~700 ha) are currently progressing towards Habitat 5330. In addition, approximately 1200 ha are occupied by active agricultural land and 600 ha by urban areas and infrastructure.

The study was conducted in the Valle della Soda area, a north–south oriented valley located in the northern sector of Monte Argentario (Tuscany, Italy). A 2.5 ha study polygon was delineated to include the valley bottom (impluvium) and adjacent slopes, encompassing the main habitat mosaic present in the area (Natura 2000 habitats 9340, 5330 and 6220*) (Fig. 2a, b); the site falls within the officially designated “demarcated area” according to EU Regulation 2020/120 (Article 4).

To determine the presence or absence of Xylella fastidiosa subsp. multiplex ST87 (Xfm ST87) in the study site, annual surveys, from 2021 to 2025, were carried out following EFSA and EPPO guidelines (EFSA 2023; EPPO 2023) and in accordance with the Commission Implementing Regulation (EU) 2020/1201. The area spans an elevational gradient from 130 to 240 m a.s.l., allowing the inclusion of structurally heterogeneous vegetation types within a limited spatial extent. An integrated field–UAV surveys approach was adopted.

In May–June 2021, a complete floristic survey of the study area was conducted to document the species pool of the site. In the same period, vegetation surveys were performed within permanent plots. Twelve permanent plots were established to sample the habitat mosaic occurring within the study area. Plot size was defined according to physiognomy: 10 × 10 m in tree-dominated stands and 5 × 5 m in shrub- and herb-dominated stands. Plots were georeferenced and permanently recorded to ensure precise relocation during the second survey campaign in 2025 (Fig. 2c ; Suppl. material 1: table S3). The same plots were resurveyed in June–July 2025 under comparable phenological conditions.

Field surveys allowed confirmation of the taxonomic identification of species recognizable on the orthomosaics and the recording of taxa not detectable from aerial imagery due to canopy overlap, small size, or chromatic similarity. Species cover within each plot was quantified through mapping on ultra-high-resolution RGB orthomosaics derived from UAV imagery, allowing precise spatial delineation of vegetation patches without relying on visual field estimates.

Botanical nomenclature follows Bartolucci et al. (2024) and subsequent updates provided by the Portal to the Flora of Italy (Martellos et al. 2020). Biological life forms and chorological types were assigned following Pignatti et al. (2017–2019). Diagnostic species for Annex I habitats (Directive 92/43/EEC) were specifically recorded to support habitat attribution and vegetation classification.

UAV surveys were conducted in both vegetation-sampling years using different UAV platforms. In 2021, flights were performed with a DJI Phantom 4 Pro equipped with a 20 MP RGB sensor (1" CMOS; 24 mm equivalent focal length). In 2025, surveys were conducted using a DJI Mini 4 Pro equipped with a 48 MP RGB sensor (1/1.3" CMOS; 24 mm equivalent focal length). Automated flight missions were programmed at a constant elevation of 25 m above ground level, with 60% forward overlap and 80% side overlap. Based on flight elevation and sensor geometry, the theoretical Ground Sampling Distance (GSD) was approximately 0.68 cm/pixel for the Phantom 4 Pro and 0.44 cm/pixel for the Mini 4 Pro under ideal geometric conditions. A summary of the UAV survey parameters for both campaigns is provided in Table 1.

No ground control points (GCPs) were deployed; georeferencing relied on the onboard GNSS system of the UAVs. Given the limited spatial extent of the study area (2.5 ha), the ultra-high spatial resolution of the imagery, and the plot-based analytical framework, relative spatial consistency between survey years was considered sufficient for accurate intra-plot comparisons.

Flights in both years were conducted under stable weather conditions, with clear or slightly hazy skies, no precipitation, and low wind speed (< 5 m/s), to minimize motion blur and illumination variability. Missions were carried out between 10:30 and 14:30 local time to ensure high solar elevation and reduced shadow effects. Camera settings were maintained consistent across flights using automatic exposure and fixed white balance to minimize chromatic variability between survey years.

Image alignment, dense point cloud generation, orthomosaic production, and 3D reconstruction were performed in Agisoft Metashape (v. 1.7.5) (Agisoft LLC 2021) (Fig. 3). Model quality was assessed using the reprojection Root Mean Square Error (RMSE), which averaged 0.29 pixels. Orthomosaics were georeferenced in WGS 84/UTM zone 32N (EPSG: 32632) and analyzed in QGIS 3.42.0.

Within each permanent plot, species cover was quantified through the manual delineation of individual taxa on ultra-high-resolution RGB orthomosaics, performed by a human operator.

Only species previously recorded and taxonomically confirmed during field surveys were mapped. Field observations conducted in both years were used to validate species presence and support image-based delineation, thereby reducing misclassification among morphologically similar taxa. Mapping focused primarily on dominant or structurally distinguishable species exhibiting consistent canopy texture or chromatic patterns. Absolute cover values (m²) were converted into percentage cover as (species area / plot area) × 100, allowing direct interannual comparison (2021–2025).

To assess the spatial distribution of plant species and validate the adequacy of the sampling effort, the Species-Area Relationship (SAR) (Tjørve and Tjørve 2017) was modelled using the Arrhenius power function (S = c A^z), where S is the number of species, A is the area, c is a constant representing the expected richness per unit area, and z represents the species accumulation rate (Arrhenius 1921).

The SAR was fitted using a log-log linear regression on 13 data points, comprising the 12 permanent plots individually (areas of 0.0025 ha for shrub- and herb-dominated stands and 0.01 ha for tree-dominated stands, with species richness recorded during the 2021 survey) and the total study polygon (2.5 ha, with the total number of species recorded through the complete floristic survey conducted in 2021). The Arrhenius power function S = c · A^z was linearized by log-log transformation for regression fitting, with the z-value estimated as the slope of the ordinary least squares regression and c derived from the intercept.

It should be noted that the plot-scale areas correspond to spatially non-contiguous permanent plots rather than nested sampling units, and that the two plot sizes reflect differences in vegetation physiognomy rather than a systematic area gradient. The resulting exponent therefore reflects scaling between individual sampling units and landscape extent rather than a classical nested-area design, and direct comparisons with z-values derived from strictly nested designs should be made with caution.

Additionally, to test the ecological plausibility of the observed richness, the expected number of species at the landscape scale was independently calculated using literature-based z-values typical for Mediterranean shrublands (Drakare et al. 2006; Tjørve 2009).

Vegetation data were organized as a species × plot matrix using percentage cover values recorded in 2021 and 2025. To evaluate potential temporal shifts in community composition, analyses were conducted on square-root transformed cover data (Legendre and Gallagher 2001). The square-root transformation was applied to reduce the disproportionate influence of highly dominant taxa, particularly Ampelodesmos mauritanicus, which structurally characterizes parts of the study area. This transformation allowed greater sensitivity to change in subdominant and shrub-layer species.

Floristic–vegetational patterns among permanent plots were explored through hierarchical cluster analysis. Species cover values were square-root transformed to reduce the influence of dominant taxa while preserving relative abundance structure. Dissimilarity between plots was calculated using the Bray–Curtis index, which is widely recommended for ecological community data due to its sensitivity to compositional differences and its robustness to joint absences (Legendre and Legendre 2012).

Agglomerative hierarchical clustering was performed using the unweighted pair-group method with arithmetic mean (UPGMA; average linkage). The resulting dendrogram was visually inspected and subsequently partitioned into alternative cluster solutions (k = 2 and k = 3) to explore primary and secondary compositional structures. Cluster boundaries were identified by examining the distribution of fusion heights in the dendrogram. All multivariate analyses were conducted using the Bray–Curtis dissimilarity matrix derived from transformed cover data.

To evaluate temporal changes in local plant biodiversity, alpha-diversity was assessed for each permanent plot by calculating Species Richness (S), the Shannon Diversity Index (H’) (Shannon 1948), Simpson’s Index (1−D) (Simpson 1949), and Pielou’s Evenness (J’) (Pielou 1966). Species richness was defined as the total number of vascular plant species recorded per plot. Shannon diversity was calculated as: H’ = −Σ p_i ln(p_i) where p_i represents the proportional cover of species i within each plot. All indices were calculated from raw percentage cover values, consistently with standard alpha-diversity practice. To account for the repeated-measures design of the permanent plots, differences in alpha-diversity metrics between the first survey (2021) and the resurvey (2025) were assessed using the Wilcoxon signed-rank test because the data were not normally distributed.

To visualize floristic gradients and assess the spatio-temporal dynamics of the plant communities, unconstrained ordination was performed using Principal Coordinates Analysis (PCoA). The ordination was based on the same Bray–Curtis dissimilarity matrix applied in the hierarchical clustering, derived from square-root transformed species cover data. PCoA was performed to visualize multivariate differences among plots and sampling years in reduced ordination space. Temporal separation was assessed by inspecting clustering patterns along the first two principal coordinate axes.

To test for statistically significant differences in species composition, a Permutational Multivariate Analysis of Variance (PERMANOVA; Anderson 2001) was conducted. For PERMANOVA, pseudo-F statistics, R² values, and p-values were reported.

To test for temporal changes in community composition, Year (2021 vs. 2025) was used as the explanatory factor. To account for the repeated-measures design of the permanent plots, permutations were constrained within each Plot ID (used as a block/stratum), ensuring that the temporal effect was tested within the same spatial units. A total of 9999 permutations was specified. A complementary PERMANOVA model including Plot as an explicit factor was also fitted to quantify spatial structuring among permanent plots. Cluster membership derived from hierarchical clustering was used for descriptive interpretation of vegetation assemblages and was not included as an explanatory factor in the PERMANOVA models. Homogeneity of multivariate dispersion among years was evaluated using PERMDISP (Anderson 2006) to ensure that significant PERMANOVA results, if detected, were attributable to shifts in centroid location rather than differences in within-group variability.

Similarity Percentage (SIMPER) analysis (Clarke 1993) was performed to identify the species contributing most to the average Bray–Curtis dissimilarity between surveys conducted in 2021 and 2025. The analysis was applied to square-root transformed species cover data, consistent with the transformation used for distance-based multivariate analyses.

SIMPER decomposes the overall between-group dissimilarity into the contribution of individual species by calculating the average contribution of each taxon to the Bray–Curtis distance between groups. Species contributions were expressed both as absolute average dissimilarity and as percentage contribution to total between-year dissimilarity. Cumulative contributions were used to identify the set of taxa accounting for the majority of interannual compositional differences.

Flight missions were planned and executed using the DJI Fly® application (DJI 2024). UAV images were processed through a photogrammetric workflow in Agisoft Metashape Professional® (64 bit) to generate dense point clouds and orthomosaics.

Spatial analyses, vegetation interpretation, and cartogra­phic outputs were performed using QGIS (version 3.42.0) (QGIS Development Team 2024).

Statistical analyses were conducted in R (R Core Team 2024), and community-level metrics were calculated using the vegan package (Oksanen et al. 2023).

The vegetation dataset supporting the results of this study is available as Supplementary material (Suppl. material 1: table S2).

During the monitoring, 34 taxa were analyzed as a whole, comprising a total of 454 samples. The presence of Xfm ST87 in the study site could be confirmed in each year, although estimated annual incidence varied considerably: the lowest was in 2022 (0.062), when the sampling effort was greatest (192 analyzed), and the highest was in 2024 (0.186), when 43 plants were analyzed. This variation in sampling efforts reflects both the need, in the early years, to test a broader range of plant species in order to identify those in which the bacterium could most easily be detected by PCR, and the specific aim of this monitoring, namely to confirm the consistent presence of the pathogen over time rather than to quantify its prevalence within the study area (Suppl. material 1: fig. S1).

A total of 90 vascular plant species were recorded within the study area, including both the permanent survey plots and the surrounding area (Suppl. material 1: table SS1).

The overall floristic composition reflects the Mediterranean sclerophyllous vegetation framework typical of Natura 2000 Habitat 5330 (Thermo-Mediterranean and pre-desert scrub), with localized structural elements referable to Habitat 9340 (Quercus ilex forests).

Habitat 6220* (Pseudo-steppe with grasses and annuals of the Thero-Brachypodietea) occurs in small patches within the broader landscape matrix but is only marginally represented within the sampled permanent plots (Suppl. material 1: table S2).

The biological spectrum highlights a strong structural dominance of woody life forms, reflecting the physiognomy of the local habitats. Phanerophytes (including nanophanerophytes) constitute the most abundant life form, accounting for 58.9% of the total flora, followed by hemicryptophytes (13.3%), chamaephytes (11.1%), and therophytes (11.1%).

From a phytogeographical perspective, Mediterranean chorotypes dominate the assemblage (63.3% of species), including 44.4% Steno-Mediterranean and 18.9% Euri-Mediterranean elements. Eurasian and paleotemperate chorotypes represent a minor fraction (13.3%), mainly linked to the shadier understory conditions of Habitat 9340 or localized mesic microhabitats.

UAV imagery integrated with GIS-based spatial analysis was employed to quantify absolute and relative ground cover and to assess diachronic changes between 2021 and 2025 (Suppl. material 1: fig. S2).

Field surveys remained essential for the accurate identification of most herbaceous species and small shrub or tree seedlings. Drone imagery proved fully reliable for larger shrub and arboreal taxa.

Within the permanent plots monitored in 2021 and 2025, a total of 43 species were identified (Suppl. material 1: table S2).

Vegetation is primarily composed of two dominant taxa showing the highest cover and frequency values: Ampelodesmos mauritanicus (mean cover 37%; 100% frequency) and Erica multiflora (mean cover 14%; 80% frequency), confirming the prevalence of shrub–grass typical of Habitat 5330.

A log-log linear regression of species richness against area was applied across the 13 data points (12 permanent plots and the total study polygon) to estimate the general species accumulation trend (R² = 0.886, p < 0.001). The estimated scaling exponent was z = 0.299 (± 0.032 SE) and the intercept parameter c = 63.49, giving the species–area function:

S = 63.49 × A^0.299

(with A expressed in hectares).

Given the non-nested structure of the design (see Methods), these statistics should be interpreted in light of their primary purpose: assessing whether the observed species richness of 90 taxa at 2.5 ha is ecologically consistent with theoretical expectations.

The observed value is broadly consistent with the model prediction (S = 83.5) and falls within the range predicted from literature-based z-values for continuous Mediterranean vegetation (z = 0.15–0.25; predicted range: 73–80 species), supporting the adequacy of the sampling design (Suppl. material 1: fig. S3).

Cluster analysis applied to vegetation survey data identified three compositional assemblages (Cluster A, Cluster B and Cluster C) within the study area (Suppl. material 1: table S2; Fig. 4). At a Bray-Curtis dissimilarity threshold of approximately 0.54, the dendrogram clearly resolves these three distinct groupings, indicating meaningful floristic differentiation within the sampled vegetation mosaic.

Cluster A (Plot 6, 7, and 8) is characterized by higher cover of woody sclerophyllous taxa including Arbutus unedo, Pistacia lentiscus, Erica arborea, and Quercus ilex, species diagnostic to Habitat 9340.

Cluster B (Plot 1–5, 9, 11, and 12) encompasses the majority of surveyed plots and is dominated by Ampelodesmos mauritanicus, Erica multiflora, Cytisus laniger, and various Cistus species, diagnostic to Habitat 5330.

Cluster C (Plot 10) captures a subgroup floristically and physiognomically attributable to an intermediate vegetation type between Habitats 5330 and 9340.

Cluster analysis consistently grouped surveys conducted within the same plot in 2021 and 2025 into paired subclusters. Each 2021 survey clustered with its corresponding 2025 survey, indicating strong compositional consistency between the two survey campaigns. Observed interannual differences were limited to minor shifts in relative dominance among co-occurring woody species, with no evidence of directional floristic change.

Species richness (S), Shannon Diversity (H’), Simpson’s Index (1−D), and Pielou’s Evenness (J’) were compared between 2021 and 2025 using Wilcoxon signed-rank paired tests to account for the repeated-measures design (n = 12 plots).

Species Richness (S) did not differ significantly between years (Wilcoxon signed-rank test, W = 5.0, p = 0.750). The Shannon Diversity Index (H’) showed a marginal numerical variation, shifting from a mean of 1.690 (± 0.327 SD) in 2021 to 1.711 (± 0.343 SD) in 2025 (W = 13.0, uncorrected p = 0.043).

Similarly, Simpson’s Index (1−D) and Pielou’s Evenness (J’) displayed positive trends just below the significance threshold (1−D: W = 15.0, p = 0.064; J’: W = 16.0, p = 0.077). As median and interquartile ranges were broadly comparable between survey years (Fig. 5), these subtle fluctuations represent microscopic adjustments in relative species abundances rather than robust structural changes, confirming the overall stability of local alpha-diversity over the study period.

The first two axes explained 41.8% and 16.4% of the total variation, respectively (58.2% cumulative). Ordination revealed strong overlap between sampling years, with each permanent plot maintaining close spatial proximity to its corresponding resurvey. No consistent directional displacement between years was evident along the main ordination axes, and paired surveys clustered tightly in multivariate space.

To formally test temporal effects while accounting for the repeated-measures design, a PERMANOVA was performed with permutations constrained within plots (strata = Plot). The analysis detected a statistically significant effect of Year (Pseudo-F = 0.087; p = 0.007). However, the explained variance associated with the temporal factor was extremely low (R² = 0.0039), indicating that Year accounted for less than 0.4% of total compositional variability (Fig. 6).

A complementary model including Plot as an explicit factor confirmed significant spatial structuring among plots (R² ≈ 0.995; p < 0.001), whereas the marginal contribution of Year remained negligible.

No significant differences in dispersion between years were detected (F = 0.0011; p = 0.972), indicating comparable within-group variability and excluding dispersion-driven effects in the PERMANOVA results.

Collectively, multivariate analyses indicate strong spatial differentiation among plots and minimal temporal compositional change between 2021 and 2025.

Although the overall interannual dissimilarity was low and the temporal effect negligible in PERMANOVA terms (R² = 0.0039), SIMPER analysis was performed to identify which species, and in particular which Xylella fastidiosa target species, contributed most to the observed between-year compositional differences. This information is directly relevant to the phytosanitary monitoring objectives of the study, as it allows assessment of whether the minimal floristic change detected was concentrated in susceptible taxa or distributed diffusely across the community.

Overall interannual dissimilarity was low, consistent with the limited variance explained by Year in PERMANOVA. The highest individual contributions to between-year dissimilarity were associated with Cytisus laniger (36.3%), followed by Brachypodium sylvaticum (11.6%), Erica multiflora (8.3%), and Cistus monspeliensis (6.5%) (Fig. 7). Notably, Cytisus laniger, the Xfm ST87 target species showing the most marked cover increase over the monitoring period, accounted alone for over one third of total interannual dissimilarity, reflecting its floristic dynamism within an otherwise stable community. The remaining target species (Spartium junceum, Rhamnus alaternus, Cistus monspeliensis) contributed marginally and showed no directional change. Between-year dissimilarity was otherwise distributed diffusely across common perennial taxa, with no single non-target species accounting for a disproportionate share of variation.

The total vascular flora recorded in the Valle della Soda study area, 90 species over a 2.5 ha polygon, reflects the floristic richness expected for well-preserved Mediterranean maquis and garrigue mosaics of the Tyrrhenian coast. The dominance of Steno-Mediterranean and Euri-Mediterranean chorotypes (63.3% combined) is consistent with the biogeographic positioning of Monte Argentario within the Mediterranean Biogeographic Region and aligns with the phytosociological framework of the Quercetea ilicis class, which characterizes the sclerophyllous evergreen vegetation of the central-western Mediterranean basin (Blasi et al. 2010; Biondi et al. 2014). The strong prevalence of phanerophytes (58.9%), with therophytes accounting for only 11.1% of the total flora, indicates a structurally mature community at a shrubland-to-woodland transitional stage, little affected by recent large-scale disturbances such as fire or intensive grazing. The very low therophyte cover within permanent plots (below 1%) reinforces this interpretation, as annual species tend to expand under post-disturbance or soil-degradation conditions in Mediterranean ecosystems (Lloret et al. 2005).

The species-area scaling exponent obtained (z = 0.299 ± 0.032 SE, R² = 0.886, p < 0.001) falls within the theoretical range proposed for continuous, non-isolated continental ecosystems (z = 0.15–0.30; Rosenzweig 1995), and is consistent with values reported for Mediterranean shrubland and semi-natural vegetation systems (Drakare et al. 2006; Tjørve 2009), albeit at the upper boundary of this interval. The high goodness of fit confirms that the Arrhenius power function adequately describes species accumulation across the sampled spatial extents, from individual plots to the full study polygon. This value is also fully coherent with the physiognomic and structural characteristics of the Habitat 5330-dominated mosaic: as demonstrated by Keeley and Fotheringham (2003), shrub-dominated communities characterized by a dense matrix of perennial woody species and a low occurrence of annuals typically exhibit stable and relatively low z-values — precisely the pattern observed here, where phanerophytes and chamaephytes account for nearly 70% of the biological spectrum while therophytes are marginally represented.

The exponent obtained is slightly higher than broader Italian baseline estimates (z ≈ 0.26) (D’Antraccoli et al. 2024). Rather than reflecting a scale-dependent dampening effect, this modest upward deviation is best interpreted as a consequence of the structural and floristic heterogeneity of the three-habitat mosaic (H9340, H5330, H6220*), which generates higher compositional turnover across short distances than would be expected in a more physiognomically homogeneous system. In such mature sclerophyllous systems, spatial heterogeneity is driven by the gradual structural maturation of dominant taxa rather than by sharp compositional discontinuities or pioneer species invasion.

The observed richness of 90 species at 2.5 ha slightly exceeds the value predicted by the fitted model (S = 83.5) and the interval obtained from literature-based z-values for continuous Mediterranean vegetation (z = 0.15–0.25; predicted range: 73–80 species) yet remains fully coherent with SAR expectations for Mediterranean shrubland landscapes. The species richness values recorded within individual permanent plots are consistent with the fitted species-area function across the full range of sampled areas (0.0025–0.01 ha), with no systematic deviation detected at plot scale. The SAR analysis therefore serves a dual interpretative role: it validates the adequacy of the sampling design, confirming that the selected plots proportionately captured the core diversity of the 2.5 ha area, and provides a spatial framework that helps contextualize the strong compositional inertia and temporal stability observed in the multivariate analyses between 2021 and 2025.

The hierarchical cluster analysis identified a clear vegetation composition within the investigated vegetation mosaic. The two-cluster solution defines the primary organizational pattern of the system, with a secondary partition emerging at k = 3 and clearly resolved at a Bray-Curtis dissimilarity threshold of approximately 0.54. Although cluster boundaries are gradual rather than sharp, as expected in Mediterranean sclerophyllous landscapes characterized by continuous structural transitions among vegetation phases (Whittaker 1967; Legendre and Legendre 2012), the three groupings are consistently supported by the PCoA ordination results.

Cluster A, dominated by sclerophyllous arboreal taxa (Arbutus unedo, Pistacia lentiscus, Erica arborea, Quercus ilex), corresponds structurally to the Habitat 9340 woodland matrix, whereas Cluster B, encompassing the majority of surveyed plots, reflects the shrub-grass maquis physiognomy of Habitat 5330, with Ampelodesmos mauritanicus and Erica multiflora as structural dominants. Cluster C captures a secondary level of structural heterogeneity, representing a floristically transitional subgroup characterized by increased dominance of woodland-associated species while retaining shrub physiognomy. However, habitat classification under Natura 2000 requires concordance between floristic composition and physiognomic structure (European Commission 2013): although woodland-associated species contribute substantially to this assemblage (~36% of mean total cover), canopy continuity remains incomplete and shrub dominance persists. The identified unit therefore represents a structurally mature expression within Habitat 5330, showing woodland affinity without reaching the structural threshold of Habitat 9340.

No independent therophyte-dominated assemblage emerged from the multivariate analysis, and no cluster was characterized by dominance of annual species. The very limited representation of therophytes (mean cover below 1%) indicates that Habitat 6220* elements, while potentially present in the broader landscape matrix, do not constitute a floristically coherent unit within the sampled area. The floristic-vegetational data collected are consistent with those reported by Arrigoni and Di Tommaso (1997) for similar environmental settings on Monte Argentario.

Temporal consistency reinforces this structural interpretation: surveys conducted in 2021 and 2025 consistently grouped within the same clusters, with each plot pairing its two survey years into tight subclusters. Observed changes were restricted to compensatory dominance shifts among co-occurring woody taxa, without evidence of directional turnover. Such dynamics are characteristic of mature Mediterranean sclerophyllous systems, where structural maturation and canopy development prevail over rapid compositional replacement (Keeley and Fotheringham 2003). Overall, the convergence of hierarchical clustering, ordination analysis, and temporal stability indicates that the study area constitutes a spatially coherent and structurally mature Mediterranean shrub–woodland mosaic, whose detected differentiation reflects internal successional gradients within Habitat 5330 rather than discrete habitat entities, consistent with ecological theory describing community structure as a continuum shaped by gradual environmental and structural transitions (Whittaker 1967; McGill et al. 2006).

The most striking result of the present study is the remarkable compositional stability of the monitored vegetation between 2021 and 2025. Although the PERMA­NOVA formally detected a significant temporal effect (p = 0.007), the associated explained variance was negligible (R² = 0.0039), indicating that Year accounted for less than 0.4% of total floristic variability. This outcome illustrates the well-documented distinction between statistical significance and ecological relevance in multivariate community analyses (Anderson 2017; Legendre and Legendre 2012): with permutations constrained within plots, even microscopic shifts in relative cover become detectable, yet their biological meaning remains trivial. The much larger contribution of Plot identity (R² = 0.995) confirms strong spatial differentiation among plots, consistent with the local environmental heterogeneity of the study area. The PCoA ordination is fully consistent with this interpretation. Paired surveys from the same plot clustered tightly in multivariate space, with trajectories so short as to be indistinguishable from sampling noise, and no consistent directional displacement was detectable along either principal axis.

SIMPER analysis corroborated this picture while adding a relevant layer of interpretation. The decomposition of the negligible interannual dissimilarity revealed that Cytisus laniger alone accounted for 36.3% of between-year floristic change, a disproportionate contribution for a single taxon in an otherwise stable community, directly reflecting its significant cover increase documented in the Target dynamics section. This result confirms that the minimal compositional shift detected between 2021 and 2025 was not distributed randomly across the flora but was concentrated in the most dynamically active Xfm ST87 target species. The remaining target species contributed marginally to interannual dissimilarity, consistent with their cover stability over the monitoring period. Together, these analyses paint a coherent picture of strong structural and compositional inertia, a hallmark of mature, stress-tolerant Mediterranean woody vegetation (Blondel and Aronson 1999; Lloret et al. 2012), while identifying Cytisus laniger as the floristic element most deserving of attention in future resurveys.

Such structural inertia is not merely a passive property of established vegetation. Mediterranean woody communities, and maquis formations in particular, are known to exhibit considerable buffering capacity against moderate disturbance through resprouting, clonal growth, and the persistent soil seed bank (Pausas et al. 2004; Verdú and Pausas 2007). The absence of any detectable floristic collapse over the monitoring period therefore does not preclude the ongoing presence of Xfm ST87 in the study area, but indicates that the plant community as a whole has not yet experienced measurable phytopathological impact at the community scale. Landscape-level effects of Xylella fastidiosa subsp. pauca on olive groves in Apulia have been documented only after a decade or more of progressive symptom development (Saponari et al. 2019; Boscia et al. 2017), suggesting that the period examined here may be insufficient to detect early-stage floristic shifts in species with slow demographic turnover such as the woody taxa dominating this study area.

Alpha-diversity results are consistent with the overall picture of community stability, while revealing subtle internal shifts in species dominance. Species richness remained perfectly stable between the two survey campaigns (W = 5.0, p = 0.750), indicating that no local extinctions or colonizations occurred within the permanent plots, a result consistent with the conserved floristic matrix documented by multivariate analyses.

While the Shannon Diversity Index displayed a slight numerical fluctuation between the two survey campaigns (from 1.690 to 1.711; uncorrected p = 0.043), alongside marginal positive trends for Simpson’s Index (1−D) and Pielou’s Evenness (J’), the absolute magnitude of these shifts was ecologically negligible (delta = +0.021 for Shannon). Thus, the overall pattern is best interpreted as stable: the community did not genuinely diversify, but the distribution of relative abundances shifted very slightly towards greater evenness.

This subtle shift in uniformity, in the absence of any change in species richness, is consistent with a partial redistribution of cover between species rather than compositional substitution. For example, the slight decline in cover dominance of Ampelodesmos mauritanicus is offset by the increase in subdominant species. Such minor fluctuations in diversity metrics without changes in species richness have been documented in Mediterranean shrublands undergoing post-disturbance recovery or interannual climatic variability (De Cáceres et al. 2015; Lloret et al. 2012), and do not necessarily imply any directional ecological trend.

The analysis of the Xfm ST87 target species recorded in the permanent plots (Cytisus laniger, Spartium junceum, Rhamnus alaternus and Cistus monspeliensis) constitutes perhaps the most ecologically significant component of the present study, as it directly addresses the primary concern motivating the monitoring design.

None of the target species showed a decline in cover over the study period. Spartium junceum and Rhamnus alaternus were perfectly stable across all twelve plots (Δ = 0 in all cases). Cistus monspeliensis was largely stable across plots (Δ = 0 in 10 out of 12 plots), with only minor and non-directional cover fluctuations (contributing 6.5% to interannual dissimilarity). Cytisus laniger, one of the Xylella fastidiosa subsp. multiplex ST87 host species reported in the Mediterranean maquis of central Italy (EFSA 2026; Falsini et al. 2025), showed not a decline but a statistically significant increase in mean cover, rising from 6.37% in 2021 to 10.97% in 2025 (W = 0, p = 0.0078), with cover increasing in eight out of twelve plots and remaining unchanged in the remaining four. This result is unambiguous and counter-intuitive from a strictly pathological perspective.

Several non-mutually exclusive interpretations can be advanced. First, the absence of visible decline in C. laniger may reflect the early and sub-symptomatic phase of infection, during which bacterial colonization progresses through the xylem vessels without yet triggering the severe water stress and scorching symptoms that characterize advanced disease stages (Saponari et al. 2013; Sicard et al. 2018). Second, it is possible that the local populations of C. laniger exhibit a degree of tolerance or partial resistance to the bacterium, as has been documented for certain genotypes of Xylella hosts in other contexts (Velasco-Amo et al. 2020). Third, the observed cover increase may simply reflect the intrinsic growth dynamics of a fast-growing, resprouting leguminous shrub, whose canopy expansion between the two survey campaigns is well within the range of natural interannual variability in Mediterranean shrublands, independent of any pathological process. Finally, the monitoring window may simply be insufficient to detect early-stage population declines in a perennial resprouting shrub, whose vegetative growth can mask underlying root and xylem damage for several growing seasons before external symptoms become visible (Saponari et al. 2019).

Regardless of the mechanism involved, the absence of vegetation variation over the 2021–2025 period is an important baseline finding. It provides a reference against which future resurveys, particularly if extended to a ten-year horizon, can detect the onset of plant landscape change.