Phytoindication as a Unique Ecological Tool for Assessing the Response of Biotic Systems to Environmental Factors

1. Introduction

Global climate change is profoundly altering key ecological processes and the conditions necessary for the survival of living organisms [1]. Traditional concepts of environmental factors, developed throughout ecological history, must be reconsidered, as the patterns of their formation and interaction under contemporary changes differ significantly from those underpinning conventional interpretations. Anthropogenic pressure has reached such a magnitude that even factors and complexes previously regarded as natural, such as climate, soil-forming rocks, topography, biota, and the progression of ecological time, have transformed [2,3]. Today, virtually no regions on the planet remain free from human influence. This situation imposes new demands on our understanding of the role of biodiversity reserves: they should not only be strictly protected areas but also include any regions capable of ensuring the conservation of rapidly vanishing biodiversity.

Ecology conceptualises abiotic factors as distinct components of the physical and chemical environment that influence organisms' distribution, abundance, and survival within an ecosystem [4,5]. The attribute “distinct” implies a conceptual interpretation of factors as theoretically independent phenomena. However, in real ecosystems, these factors often interact or are interrelated, and complete independence is seldom observed. Therefore, in practical research, environmental gradients are used as proxies for environmental factors. An environmental gradient represents an environmental factor as an ordered value across space or time, reflecting the degree of differences between locations or moments in time along the continuum of this factor [6,7]. This idea was reflected in the concept of the vegetation continuum represented along an environmental gradient, which formed the basis for gradient analysis [8]. In turn, gradient analysis of vegetation is a development of the ideas of the factorial approach for explaining soil differentiation [9,10,11]. Soil properties are a function of five environmental factors: parent material, regional climate, topography, biota, and time [12]. These five soil-forming factors are clearly the same factors that determine vegetation. Vegetation distribution corresponds to regional climatic patterns, varies with changes in soil parent material, differs on slopes with opposite aspects, changes according to geographic differences between floras, and evolves [13]. Environmental gradients are classified into indirect (complex), direct, and resource gradients [14]. Indirect gradients, such as elevation above sea level, influence plant growth indirectly by correlating with other environmental factors like precipitation, wind, and temperature. Direct gradients exert a direct physiological effect on plants, while resource gradients consist of factors organisms directly utilise as resources for growth. There is no clear boundary between direct factors and resources. Under certain conditions, a resource may function as a direct factor [15,16]. Plant growth performance can be represented as a function of four main environmental factors: nutrients, water, temperature, and light [10,17], which correspond to the primary indicator values in virtually all indicator systems.

Plants can serve as indicators of environmental factors [18,19]. The conceptual foundation for applying bioindication is the theory of species' ecological niches [20]. This theory posits that the species composition of a given community allows for the assessment of environmental conditions [21]. The bioindication process relies on indicator values that represent the ecological optimum of a species along the gradient of a specific environmental factor [22]. These values are incorporated into a community-weighted mean calculation based on cover-abundance scores for the relevant plant community [23,24]. This calculation is then used to indicate the level of the corresponding environmental factor under which the community has developed [25]. Phytoindication enables the assessment of various physical abiotic environmental factors, including site illumination, temperature regime, climate continentality, soil or substrate moisture levels, acidity (pH), nutrient availability, and, in some cases, salinity [26,27]. Indicator scales can be used to interpret the results of ordination analysis, as they allow the identified gradients in species composition to be linked to the abiotic environmental factors that influence them [28]. Gradient analysis and biological indication share a common theoretical foundation and can be integrated effectively. Weighted average indicator values calculated for sample plots can be overlaid onto the ordination space as passive variables or used to construct isolines representing the direction and intensity of environmental gradients [29,30]. Phytoindication estimates of environmental properties serve as auxiliary variables in ordination analyses, facilitating the ecological interpretation of the primary gradients shaping plant community structure [31]. This approach enables quantification of the positions of studied taxa along moisture gradients, soil pH, trophic status, and light availability without directly measuring abiotic factors. It allows for statistical testing of differences between plant communities [26].

Phytoindication integrates environmental conditions over extended periods, eliminating the need for repeated and costly measurements, and can be applied in situations where direct assessments are impossible [24,32]. It enables the retrospective reconstruction of environmental parameters based on species composition, while species response models allow for the calibration of abiotic indicators [25,33,34]. In addition to natural factors, levels of anthropogenic transformation are indicated using scales of naturalness and hemeroby [35,36]. Standardised systems, such as the Ellenberg scales, ensure comparability of results across regions and studies by relying on data from standardised geobotanical relevés [37]. When combined with geostatistical methods and remote sensing, indicator values facilitate the creation of spatial models of environmental factors, which is especially valuable in large-scale studies where direct measurements at each site are technically or financially unfeasible [38]. The capacity for spatial modelling of indicator values based on remote sensing data is predicated on the premise that the composition of plant species, which determines the weighted mean indicator value for a given site, is associated with the spectral characteristics of the vegetation cover [39]. The biomass structure, leaf area, chlorophyll and water content, seasonal dynamics and spatial heterogeneity of different species assemblages are all subject to variation, shaped by specific abiotic conditions [40]. These variations can be reflected in spectral indices derived from satellite or aerial imagery [41]. When these spectral indicators and other spatial predictors, such as elevation, distance to water, and climatic data, are statistically linked to known field estimates of indicator values, it becomes possible to develop a model that predicts these values across an entire region [42,43]. Remote sensing facilitates the extrapolation of bioindication data from specific points to broader areas, thereby ensuring the acquisition of high-resolution spatial information [44].

A fundamental challenge in interpreting which indicator values represent specific environmental factors arises because these values do not directly measure them. Instead, they reflect the responses of individual plant species and plant communities to those factors [45,46]. Biotic systems can integrate external influences [47]. The complex filtering mechanism underpinning the formation of diverse plant assemblages in nature is a fundamental factor in the subsequent variation in indicator scores of species both within and between plots [27]. The law of the limiting factor [48,49] posits that a biotic system responds primarily to factors that are either deficient or in excess, exhibiting little or no response when a factor is within the species' optimal range. This explains why nominal indicator scales often show cross-sensitivity to multiple factors. For instance, the indicator scale for available nitrogen in the Ellenberg system has been demonstrated to correlate with the productivity of plant communities [50], as it reflects not only soil nitrogen content but also other fertility factors such as phosphorus [22,51]. Consequently, in recent adaptations of the Ellenberg indicator system, this scale has been renamed the nutrient availability scale [22]. Similarly, the indicator values for soil reaction reflect changes in pH and the quantity or saturation of exchangeable Ca²⁺ [52]. Amplitude indicators of water regime dynamics show a stronger correlation with soil moisture indicator values than average water regime indicators [52]. The accuracy of these indicator scores depends on the contrast of environmental conditions, often referred to as the length of the gradient. In cases of short gradients, indicator values may be influenced more by random fluctuations in species composition than by the actual environmental gradient [25]. A significant conceptual inconsistency exists within the framework of phytoindication. This inconsistency pertains to the mechanistic transfer of autecological properties of individual plant species, as indicated by indicators of their optimal zones along environmental gradient factors, to the community level for environmental factor assessment [53]. Consequently, estimates of a plant species’ ecological niche are used as scales to indicate environmental properties at the community level. The presence of a species in a description can result from numerous factors, including environmental conditions, interspecific interactions, or neutral influences [54]. Nonetheless, it is generally accepted that ecological indicator values for a site represent an integrated signal of species-environment relationships at the community level and provide reliable information about the long-term environmental conditions characterising that site [55]. Comparisons with empirical measurements of environmental factors reveal shifts in indicator scores across different phytosociological vegetation classes [56], which may also reflect significant changes in the informational value of plant species depending on the context within different communities.

Thus, the history of phytoindication development leaves several unresolved questions: Is phytoindication a convenient and relatively inexpensive surrogate for measuring environmental factors, or is it the only viable method to assess the impact of environmental factors on plant communities? To what extent is it methodologically justified to transform factors influencing plant growth performance into primary indicator values within indicator systems? What exactly is subject to phytoindication, and to what degree can it be applied to assess ecosystem specificity and facilitate their comparison? In this study, we aimed to evaluate the response of vegetation cover to the combined influence of natural and anthropogenic factors within a protected area located in the heart of an industrial region and subject to substantial pressure on biodiversity. By applying ecological indicator values alongside assessments of naturalness and hemeroby, we examined how plant communities integrate the effects of multiple drivers. Our results demonstrate that the use of hemeroby evaluation, traditionally applied in urban or ruderal ecosystems, is also appropriate and informative for protected areas, where biodiversity is increasingly shaped by pervasive human impact. This approach provides new insights into the extent to which protected landscapes, despite their conservation status, are embedded within broader socio-ecological contexts and therefore cannot be considered isolated from regional anthropogenic pressures.

2. Materials and Methods

The study was conducted in 2024 within the terrestrial ecosystems of the Dnipro-Orilskiy Nature Reserve (Dnipropetrovska Oblast, Ukraine). The Dnipro-Orilskiy Nature Reserve was established in 1990 and encompasses an area of 3,766 hectares [57]. The reserve is situated at the centre of the industrial agglomeration formed by the cities of Dnipro and Kamianske; therefore, despite its protected status, its ecosystems are subject to considerable anthropogenic pressure (Figure 1). Furthermore, the relatively short history of the reserve precludes expectations of substantial natural recovery. Nevertheless, the reserve serves as a biodiversity hotspot within this industrial region [58,59], where two groups of ecosystems are protected: floodplain water bodies of the Dnipro River and psammophytic communities of the first above-floodplain (sandy) terrace of the Dnipro valley. This study focuses on the vegetation of the sandy terrace. A contour line at 78 metres above sea level proved to be a reliable marker of the boundary of the sandy terrace, which we used to delineate its configuration. The highest point of the terrain reaches 94 metres above sea level. The area is characterised by sandy dunes redeposited under the influence of aeolian processes. The soils are classified as Arenosols. In the northeast of the study area, the terrain is intersected by the valley of the small Protich River, where the soils are classified as Fluvisols and Cambisols.

During the summer of 2024, we recorded the presence of all vascular plant species in 1079 sample plots of 4 × 4 m. The average distance between the sites was 65 ± 34 m (in 95% of cases, it ranged from 21 to 112 metres). The area of the convex hull around the point cloud was 12,090,087 m², so the area surrounding each site averaged 11,204 m².

The results of our records are available through the GBIF database [60]. The projected species coverage was expressed as a percentage. For this study, infraspecific taxa were generally treated as species. Critical specimens were collected and identified using microscopy. The sample plots were distributed throughout the study area, and their exact locations were determined using GPS (Garmin eTrex, ±5 m). We aimed to select samples that represented the full range of community types and plot sizes within the піщанoї тераси. Plant taxonomy was based on the Euro+Med PlantBase (http://ww2.bgbm.org/EuroPlusMed). Vegetation classification was performed using TWINSPAN software [61]. The fidelity of diagnostic species for clusters was determined using the fidelity coefficient (phi coefficient), with a fidelity threshold of 25 (and 50 for highly diagnostic species). Species with a frequency of occurrence greater than 25% (or greater than 50% for highly constant species) were considered constant. In comparison, species with projected coverage exceeding 10% were dominant [62]. The phi coefficient was calculated using the indicspecies library [63]. Syntaxes were identified up to the association level; names of syntaxes are given according to the European Vegetation of Europe [64] and the Prodromus of Vegetation of Ukraine [65]. The classification of vegetation and the Prodromus are not the subject of this study; therefore, for detailed information, we refer to the data published in the Mendeley Data repository of scientific data and materials [66].

Ecological indicator values for Europe, namely light availability, temperature, soil moisture, soil reaction, and nutrient availability, were used to represent environmental factors [28]. Continentality was assessed using the original Ellenberg scales [67,68]. The Frank and Klotz scales [69] were employed to evaluate hemeroby. The original scales were converted by calculating the mean of the minimum and maximum values for each species, which were then transformed onto a 100-point scale. The weighted average of the hemeroby scores, accounting for projected plant cover, was used to characterise the hemeroby of each sample [70]. Types of plant social behaviour are based on plant species' roles within communities. They reflect how plants are connected to their habitats and the informativeness and naturalness of these connections. The properties of the types present in a community can be used to infer the richness of ecological information within the community, its stability and naturalness, the degree of niche occupation, the regeneration capacity of the community, and the extent of disturbance, transformation, or deviation from the natural state [71].

A one-way multivariate analysis of variance (MANOVA) was performed to assess the overall effect of vegetation classification on the set of ecological indicator values, with plant association as the independent factor and six Ellenberg indicator values (Light, Temperature, Continentality, Moisture, Reaction, and Nutrients) as dependent variables. The Pillai’s trace test statistic was used to evaluate multivariate significance, as it is considered the most robust and conservative among the commonly used MANOVA test statistics, particularly in cases of unequal group sizes or violations of multivariate normality. Multivariate analysis of variance (MANOVA) using Pillai’s trace was performed in R [72] with the base function manova() from the stats package.

A Detrended Correspondence Analysis (DCA) was first conducted on the species composition matrix using the decorana() function from the vegan package in R to determine the appropriate ordination method for analysing species–environment relationships. The length of the first DCA axis was 8.64 standard deviation units, indicating a long ecological gradient and a unimodal species response to environmental variation. According to the commonly accepted criteria [73,74], Canonical Correspondence Analysis (CCA) is the most appropriate ordination method if the gradient length exceeds three standard deviation units. Therefore, CCA was applied to explore the influence of environmental variables on species composition. A partial Canonical Correspondence Analysis (partial CCA) was employed to evaluate the unique contribution of separate predictors to the inertia in species composition. This approach allows for assessing the effect of a specific explanatory variable while statistically controlling for the' influence of covariates. In each model, the species composition matrix was used as the response variable, and a single variable of interest (e.g., association type, or an individual Ellenberg indicator, naturalness, and hemeroby) was specified as the constrained (explanatory) variable. The conditioning matrix using the Condition() function included the remaining set of predictors. This procedure isolates the effect of the focal variable by removing the portion of variation in species composition attributable to the controlled ecological gradients. The analysis was performed using the cca() function from the vegan package in R. Statistical significance of the models was evaluated by permutation tests with 999 permutations under the reduced model. The adjusted R² values were calculated via the RsquareAdj() function, providing a penalised estimate of the proportion of inertia explained uniquely by each explanatory variable, accounting for the complexity of the model.

Information on canopy closure of the tree layer and the percentage cover of the shrub and herb layers was subjected to principal component analysis (PCA) to avoid multicollinearity among predictors in subsequent correlation analyses. The PCA yielded three new variables whose scores were orthogonal (i.e., independent) by definition. Based on the loading patterns on the observed variables, these components were interpreted as follows: PC1 representing Canopy and Shrub versus Herb layers, PC2 representing Shrub layer dominance, and PC3 representing coherent variation in Canopy and Herb layers. This analysis aimed to eliminate multicollinearity rather than reduce the dimensionality of the predictor space; therefore, all three principal components were retained for subsequent analyses.

For spatial modelling of the variability of phytoindication estimates of ecological factors, their regression dependence on 29 Landsat spectral indices was established using the Random Forest procedure [75]. A detailed overview of the spectral index and its calculation methods is presented in our previous article [35] and further discussed in detail in the scientific protocols repository Protocols.io [76].

3. Results

A total of 297 plant species were identified within the study area. The recorded plant communities represent 24 plant associations belonging to 14 vegetation classes. A Canonical Correspondence Analysis (CCA) revealed the effect of association type, ecological factors, naturalness, and hemeroby on species composition. The global permutation test showed these predictors' highly significant impact on species composition (χ² = 6.2, F = 13.6, P < 0.001). The adjusted R², which accounts for the number of predictors and sample size, was 0.26, indicating a moderate proportion of explainable variation in plant species composition. A partial Canonical Correspondence Analysis (CCA) was performed to isolate the unique effect of associations on species composition, while controlling for Ellenberg indicator values, naturalness, and hemeroby. The model showed a significant independent effect of association type (χ² = 3.16, F = 9.37, P < 0.001). After accounting for the other covariates, the adjusted R² was 0.132, indicating that approximately 13.2% of the explainable inertia in species composition could be attributed to plant association alone. The environmental factors represented by Ellenberg indicator values, naturalness, hemeroby, and association types included as conditional variables, explained 5.2% of the total inertia in vegetation composition (adjusted R² = 0.0515, F = 10.2, P < 0.001). A one-way MANOVA was conducted to examine the effect of plant associations on a set of Ellenberg indicator values (light, temperature, continentality, moisture, reaction, and nutrient availability). The results revealed a statistically significant multivariate effect of association on the combined set of environmental indicators (Pillai’s trace = 2.7, F(138, 6330) = 36.7, P < 0.001). This reveals that the association is connected with significant differences in the overall ecological conditions inferred from the Ellenberg indicator system. The high value of Pillai’s trace (2.7) suggests a large multivariate effect size, meaning the variation in indicator values among associations is statistically significant and ecologically meaningful.

A partial CCA model, considering plant association type, naturalness, and hemeroby as constrained variables, revealed that Ellenberg indicator values could explain 3.9% of the plant community inertia (adjusted R² = 0.039, F = 10.3, P < 0.001). The influence of environmental factors was directed primarily along the first canonical axis (Figure 2). CCA1 was positively correlated with soil moisture and nutrient availability, and negatively correlated with light availability, temperature, soil reaction, and continentality. The environmental factors ' principal component analysis (PCA) extracted two principal components with eigenvalues greater than one (Table 1). Principal component 1 closely corresponded to the direction of CCA axis 1, while principal component 2 aligned with CCA axis 2. Principal component 1 represents a gradient characterised by increasing soil moisture and nutrient availability and decreasing light availability, temperature, continentality, and soil reaction (i.e., increasing acidity). Principal component 1 after accounting for the effect of association types included as conditional variables explained 5.2% of the total inertia in vegetation composition (adjusted R² = 0.052, F = 8.0, P < 0.001).

Principal component 2 captures a gradient of increasing naturalness associated with decreasing hemeroby. Principal component 2, after accounting for the effect of association types included as conditional variables, explained 5.3% of the total inertia in vegetation composition (adjusted R² = 0.053, F = 8.2, P < 0.001).

Light availability was able to explain 4.0% of the inertia in the plant species matrix (F = 46.5, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of light availability explained 0.56% of the total inertia in species composition (F = 8.3, P < 0.001). This represented 14.4% of the total inertia explained by the set of ecological predictors derived from the Ellenberg system, highlighting the relative importance of light availability in shaping plant community composition. The phytoindication estimates of light availability ranged from 5.1 to 8.5 in 95% of cases. Community affiliation based on plant association accounted for 94.1% of the variation in light availability (F = 753.9, P < 0.001). Psammophytic communities exhibited the highest levels of light availability, with phytoindication scores consistently exceeding 7.4 (Table 1). Forest communities had light availability values below 6.0. Meadow communities and shrublands occupied an intermediate position in terms of light availability.

Temperature was found to explain 2.7% of the inertia in the plant species matrix (F = 31.4, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of temperature explained 0.68% of the total inertia in species composition (F = 9.9, P < 0.001). This represented 17.4% of the total inertia explained by the ecological predictors derived from the Ellenberg system. The phytoindication estimates of the temperature regime ranged from 4.2 to 6.1 in 95% of cases. The association membership of a plant community accounted for 66.2% of the variation in the temperature regime (F = 93.1, P < 0.001). Meadow communities exhibited the lowest temperature regime values, forest communities had intermediate values, and the highest values were observed in psammophytic communities.

Continentality was able to explain 3.8% of the inertia in the plant species matrix (F = 43.1, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of continentality explained 0.78% of the total inertia in species composition (F = 11.5, P < 0.001). This represents 20.0% of the total inertia explained by the ecological predictors derived from the Ellenberg system. Continentality scores ranged from 7.0 to 13.4 in 95% of cases. The association membership of a plant community explained 86.5% of the variation in continentality (F = 302.0, P < 0.001). The association Limonio meyeri–Festucetum pseudodalmaticae and other associations occurring on elevated terrace sites, including both psammophytic communities and the shrub association Aceri tatarici–Quercetum roboris, were characterised by continentality values generally exceeding 10. Meadow and forest communities typically exhibited continentality values below 10.

Soil moisture was able to explain 3.8% of the inertia in the plant species matrix (F = 44.1, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of soil moisture explained 0.71% of the total inertia in species composition (F = 10.5, P < 0.001). This represents 18.2% of the total inertia explained by the ecological predictors derived from the Ellenberg system. Phytoindication estimates of soil moisture ranged from 1.8 to 5.5 in 95% of cases. The membership of a community in a plant association explained 89.9% of the variation in soil moisture (F = 417.5, P < 0.001). Meadow communities formed a consistent sequence ordered by increasing soil moisture: Limonio meyeri–Festucetum pseudodalmaticae → Poёtum pratensis → Junco gerardii–Agrostietum stoloniferae → Caricetum gracilis. Forest ecosystems had intermediate soil moisture values ranging from 4.5 to 5.2 (Figure 3). Communities on elevated sandy dune areas formed a decreasing sequence of soil moisture from Salici rosmarinifoliae–Holoschoenetum vulgaris, usually occurring in interdune depressions, to those with the lowest soil moisture, such as Secali sylvestri–Caricetum colchicae and Secaletum sylvestre, typically found on the very tops of sandy dunes (Figure 4).

Soil reaction explained 1.9% of the inertia in the plant species matrix (F = 21.8, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of soil reaction explained 0.66% of the total inertia in species composition (F = 9.6, P < 0.001). This represents 16.9% of the total inertia explained by the ecological predictors derived from the Ellenberg system. The phytoindication estimates of soil reaction ranged from 5.8 to 7.4 in 95% of cases. The community membership in a plant association explained 49.5% of the variation in soil reaction (F = 46.9, P < 0.001). The highest soil reaction values were characteristic of associations such as Populetum nigro-albae, Artemisio dniproicae–Salicetum acutifoliae, and Secaletum sylvestre. Although these associations differ significantly in floristic composition (Figure 3), they share a common occurrence on elevated landforms. Conversely, the lowest soil reaction values were typical of associations such as Ficario–Ulmetum minoris, Melico nutantis–Quercetum roboris, Symphyto officinalis–Anthriscetum sylvestris, and Salici rosmarinifoliae–Holoschoenetum vulgaris, which typically occur in low-lying areas of the landscape.

Nutrient availability was able to explain 3.9% of the inertia in the plant species matrix (F = 44.9, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of nutrient availability explained 0.75% of the total inertia in species composition (F = 11.0, P < 0.001). This accounts for 19.2% of the total inertia explained by ecological predictors derived from the Ellenberg system. The phytoindication estimates of nutrient availability ranged from 1.9 to 7.7 in 95% of cases. Association membership accounted for 90.2% of the variation in nutrient availability (F = 434.5, P < 0.001). Forest communities were found under relatively high nutrient availability conditions, typically exceeding 6. Psammophytic communities occupied areas with low nutrient availability, generally below 4.5. Meadow communities occupied an intermediate position.

Naturalness was able to explain 2.2% of the inertia in the plant species matrix (F = 23.0, P < 0.001). After accounting for the influence of association types and other environmental factors included as conditional variables, the pure effect of naturalness explained 0.51% of the total inertia in species composition (F = 8.3, P < 0.001). The naturalness scores of plant communities ranged from 0.5 to 5.3 in 95% of cases. Association membership accounted for 51.6% of the variation in naturalness (F = 51.3, P < 0.001). Psammophytic communities generally exhibited the highest levels of naturalness.

Hemeroby was able to explain 2.7% of the inertia in the plant species matrix (F = 30.7, P < 0.001). The pure effect of hemeroby, after accounting for the influence of association types and other environmental factors included as conditional variables, explained 0.56% of the total inertia in species composition (F = 8.3, P < 0.001). The hemeroby scores of plant communities ranged from 14.0 to 56.5 in 95% of cases. Association membership accounted for 67.6% of the variation in hemeroby (F = 98.8, P < 0.001). The highest hemeroby levels were observed in associations such as Salicetum albae and Melico transsilvanicae–Agropyretum.

The species richness of plant communities showed significant correlations with phytoindication estimates of environmental factors (Table 2). Species richness was positively correlated with light availability, continentality, soil pH, and hemeroby. Conversely, an increase in species richness was associated with decreased soil moisture, nutrient availability, and naturalness.

The relationships among the components of the horizontal structure of plant communities also correlated with phytoindication estimates of environmental factors. Increases in the canopy and shrub layers were primarily linked to higher soil moisture and nutrient availability. In contrast, an increase in herb layer density was mainly associated with greater light availability and continentality. Shrub layer dominance increased with higher hemeroby and lower naturalness of plant communities, increased light availability and decreased temperature and continentality. Additionally, coordinated variation in the canopy and herb layers was sensitive to higher hemeroby and lower naturalness of plant communities.

The correlations between the pure effects of environmental factors and the characteristics of plant communities were more selective. The pure effects of temperature, soil reaction, and naturalness exhibited a statistically significant negative correlation with plant community species richness. The ratio of the canopy and shrub layers to the herb layer was positively correlated with the pure effect of naturalness. The pure effects of light availability and continentality were significantly positively correlated only with shrub layer dominance. In contrast, the pure effects of soil reaction and hemeroby were negatively correlated with this aspect of the horizontal structure of plant communities. Additionally, the pure impact of soil moisture was positively correlated with coordinated variation in the canopy and herb layers.

The spatial variability patterns of environmental factors distinctly divided the study area into two zones: the northern and northeastern parts, and the southern and eastern parts (Figure 5 and Figure 6). Forests and meadows predominantly characterize the northern and northeastern zone, while the south and eastern zone is dominated by psammophytic grasslands or shrublands. The spatial patterns of the various factors were quite similar, showing a clear correspondence to these zones. The pure effects of environmental factors were unique, resulting in spatial variability patterns of pure effects that were largely distinct. Furthermore, the complete and pure impact of the factors differed substantially in the scale of their spatial patterns. The full effects generally manifested at a large scale, clearly dividing the territory into two markedly different regions, whereas the pure effects were more sensitive to fine-scale spatial variations.

4. Discussion

Environmental factors influencing the biotic components of ecological systems are discussed at the earliest stages of studying ecology as a science and appear intuitively understandable, as this understanding is rooted in our everyday experience. However, ecological factors become far less straightforward and obvious when considered in the context of the history of ecology. One of the key generalisations of ecology was formulated by Liebig as early as the 1840s, well before Haeckel introduced the term ‘ecology’ in 1866. Shelford’s law of tolerance introduced non-linearity into the response of biotic systems to environmental factors [77]. In their implicit form, these laws assign particular importance to context when explaining the influence of environmental factors on biotic systems. The effect of one factor depends on the quantitative parameters of other factors, which should be understood as the response of a specific species to their combined action [78].

The bell-shaped model of a species’ response to an environmental factor [79,80] suggests that the species’ sensitivity to the factor also varies with its intensity: it is minimal both at the extremes and at the optimum. The highest sensitivity occurs in the zone corresponding to the steepest slope of the species’ bell-shaped response curve. However, the effect can be qualitatively positive (on the left branch of the curve) or negative (on the right branch). The complexity of response patterns increases when the asymmetry of the response curve and its polymodality are considered [81,82], particularly as a consequence of competitive interactions among species.

This complexity indicates that “simple” measurements of environmental factors by physical or chemical methods provide only an initial pool of information. In this context, physical and chemical processes of measuring environmental properties can compete with one another, whereas their biological equivalents cannot, in principle, match them in terms of accuracy [78]. From an ecological perspective, the significance lies not so much in the physical value of an environmental factor per se, but in the response of the living system to its influence [83]. However, the undeniable advantage of the bioindication method lies in its ability to assess the response of the ecosystem to various environmental factors [84,85]. However, there is some diversity in how phytoindication scales are defined and what they precisely indicate. These scales may be described as ‘primary environmental traits’ [21], ‘surrogates for measured environmental variables’ [25], ‘abiotic environmental variables’ [22], ‘indirect assessment of site conditions’ [86], ‘local environmental conditions’ or ‘site conditions’ in term of ecological variables [28] or ‘dynamic factors’ as opposed to static factors that can be easily measured instrumentally [87]. Indicator values are regarded as estimates of the position of species optima along gradients of environmental factors [22,28]. Therefore, by definition, phytoindication aims to identify the characteristics of species’ responses to environmental factors, rather than to determine the factors themselves. In this context, bioindication should be regarded as a leading methodological approach for measuring environmental factors, precisely the response of biotic systems to physical environmental drivers. This view fundamentally differs from considering bioindication merely as a cheap and rapid proxy for estimating environmental factors. Within an ecological framework, there is no alternative to bioindication. At the same time, physical and chemical methods serve a supplementary function as external correlates for the meaningful interpretation of observed patterns. The intuitive reluctance to accept that “precise” measurements of environmental properties using chemical or physical instruments may be imprecise from an ecological point of view creates a misleading aura of optionality around bioindication methods. Of particular note is the emphasis on evaluating the response of the biotic system within the context of hemeroby, where the difficulty or even the practical impossibility of directly estimating disturbance and human impact is considered a fundamental premise [87].

The application of bioindication necessitates a reconsideration of the nature of environmental factors' effects. The effect of “primary” factors, which align closely with our intuitive understanding of ecological drivers, can rarely be isolated and, in this respect, acts as a combination of factors. The correlation between indicator values was established for individual species and integrated assessments of environmental conditions [22]. Our results indicate that the first principal component reflects a coordinated increase in soil moisture and nutrient availability, contrasted with an increase in light availability and temperature. This pattern clearly corresponds to the differentiation of the study area into floodplain sites with lower relief, which favour the accumulation of moisture and nutrients. Conversely, elevated landforms represented by sandy hills create moisture-deficient conditions, as sandy soils have low water-holding capacity but high permeability [88], leading to rapid water loss from the root zone. The high aeration of sandy soils also enhances the mineralisation of soil organic matter [89], preventing nutrient accumulation in sandy substrates. Predictably, sparse herbaceous cover is associated with high light availability [90]. Carbonates accumulate [81] in the relief depressions, whereas sandy soils are entirely devoid of them [91,92]. This also explains the coordinated dynamics of soil pH with other phytoindication indices.

The coordinated dynamics of continentality alongside other phytoindication estimates require explanation. By definition, continentality is an ecological phenomenon operating on a much larger scale than the spatial variability of different environmental factors, which are usually site-specific [32]. Continentality of climate manifests as variation in interannual temperature and precipitation regimes along the west-east gradient across Eurasia [93]. The spatial distribution of their species determines this ecological characteristic of the plant species as a whole [94]. The interpretation of continentality values may be ambiguous at the local scale due to potential correlations with various factors, including seasonal variations in temperature and precipitation, diurnal temperature fluctuations, annual minimum temperatures, and drought conditions [22]. Thus, variability in continentality within a landscape is linked to the formation of species assemblages that share a common geographical origin and, consequently, differ in origin from other assemblages. Communities on sandy hills are represented predominantly by species of Eurasian origin, which are adapted to a higher level of climatic continentality. In contrast, floodplain ecosystems comprise species with European or Mediterranean distributions, adapted to lower levels of continentality. These features have certain ecological parallels. Existence on sandy hills with sparse herbaceous cover, caused by moisture deficit and nutrient limitation, is accompanied by substantial climatic fluctuations within a single day and over the growing season. Therefore, species preadapted to higher levels of continentality gain a competitive advantage under such conditions.

It should be noted that the identified composite environmental gradient also correlates with the level of hemeroby and the naturalness of plant communities. This relationship is likely to have geographical causes. Spatially, the floodplain ecosystems are adjacent to the reserve's ruderal areas. Thus, the floodplain ecosystems function as a buffer zone, protecting the psammophytic steppe at the reserve's centre. Consequently, a pattern emerges whereby the more humid and nutrient-rich ecosystems are somewhat more hemerobic, whereas ecosystems with opposite characteristics tend to be more natural. In this context, the first principal component represents a composite ecological gradient reflecting the coordinated spatial and temporal dynamics of a set of factors perceived as “primary.” The collection of factors defining the first principal component allows us to interpret it as a gradient prioritising natural drivers influencing vegetation cover. By contrast, the second principal component highlights the predominance of factors associated with the anthropogenic transformation of plant communities. However, it should be emphasised that anthropogenic transformation is not “refined” but represents a correlated complex of variability, including natural factors. Increasing hemeroby is associated with higher light availability, which can be explained by the reduced density of the herbaceous layer in response to anthropogenic transformation. Conversely, greater soil water availability can be regarded as a factor supporting the naturalness of plant communities. In other words, it can be assumed that plant communities in arid conditions are more sensitive to anthropogenic impacts.

Inter-correlations between indicator values are discussed as a potential source of bias in the results of bioindication [21,95]. We regard this phenomenon not as an artefact that undermines the method's accuracy, but as a regular consequence of the integrated response of biosystems to the combined effects of multiple individual environmental factors. The dimensionality of the space of environmental factors is considerably greater than the potential repertoire of biosystem responses. Consequently, the dimensionality of the reaction space is reduced, which explains the observed correlations between the indicator values. Phytoindication assessments of environmental factors reveal two principal trends of variability, which inherently represent synthetic (composite) ecological factors that are challenging to interpret in terms of any single “primary” factor. Notably, the composition of primary factors constituting a synthetic factor varies regionally and depends on the landscape, thereby being unique to each specific territory. Clearly, the composition of synthetic factors will differ under differing landscape conditions. The predominant trend of environmental variability corresponds to the significant gradient of natural conditions within the study area, accounting for 67% of the variation in ecological regimes (principal component 1). The secondary trend of variability reflects the impact of anthropogenic transformation of the vegetation cover, accounting for 15% of the variation in plant communities. This ratio underscores the respective influences of natural and anthropogenic variability in shaping the vegetation cover conditions in the study area. The variability of habitat conditions (assessed via principal component analysis of phytoindication estimates) influences the variability of vegetation structure (evaluated by the contribution of these principal components to explaining community structure, after controlling for the effect of association types included as conditional variables). This comparison yields two significant findings. First, the coordinated dynamics of environmental factors (82%) explained only 10.5% of the inertia of community composition. Second, the contribution of the minor component (hemeroby/naturalness gradient) was greater in explaining the variability of vegetation structure than that of natural factors (5.3% versus 5.2%, respectively). It is evident that the correspondence of plant communities to natural conditions has developed over an extended period of vegetation history and is reflected in the structure of vegetation types. Indeed, the explanatory power of plant associations (13.2%) exceeded that of environmental factors. It should be emphasised that plant associations and higher-rank syntaxa differ in the specific ecological conditions to which they are adapted. Thus, the contribution of natural environmental factors to vegetation variability is closely correlated with plant associations. A plant association can be understood within two principal paradigms of ecology. On the one hand, according to the structuralist paradigm, it results from interspecific interactions among plants sharing similar ecological preferences. During succession, the environment transforms, causing plants within an association to exhibit comparable ecological preferences: «The vegetation is both a cause and a product of succession; it alters the habitat and the altered habitat, in turn, determines the vegetation» [96]. On the other hand, the continuationism paradigm interprets a plant association as the outcome of the spatial co-occurrence of ecologically similar species along environmental gradients: «Every species of plant is a law unto itself, the distribution of which in space depends upon its individual peculiarities of migration and environmental requirement» [97]. Both paradigms assume the ecological similarity of species within a plant association concerning their preferences for environmental conditions. Therefore, a plant association can capture a substantial proportion of the explanatory power of phytoindication assessments of environmental factors.

The gradients of naturalness and hemeroby are less specific in relation to plant associations, and thus associations account for less of their explanatory capacity. Compared to natural drivers, anthropogenic factors operate over a much shorter timescale [98]. Consequently, plant communities in a regionally natural environment have not developed response patterns specific to the level of plant associations. Moreover, sources of anthropogenic influence are not formally correlated with the spatial distribution of plant associations [99], meaning that any association may be subject to such impacts. It should be noted that the structuralist paradigm is more effective in explaining this pattern, since anthropogenic disturbance can be regarded as a factor that returns a plant community to earlier stages of successional dynamics, which, by definition, are less specific. In turn, hemeroby can be interpreted as a gradient of anthropogenic influence on plant communities and increasing disturbance that drives communities back towards earlier and less specific stages of successional dynamics.

But what is the role of the “primary factors,” which are classically well known and intuitively clear to us, and considered universal for comparing ecosystems and landscape ecological structures? The synthetic nature of prevailing environmental gradients cannot be a basis for comparing different landscape systems due to differences in their comparative frameworks. Synthetic factors may be viewed as tools for identifying the specificity of regional conditions. Still, they cannot be used to compare different landscape systems, since only comparable entities can be compared. “Primary” factors, as independent drivers, exist as environmental conditions, but their independence disappears when we consider the consequences of their action. This disappearance of the specific effects of ecological factors is a consequence of Shelford’s law of tolerance. When species occur within their ecological optimum, their sensitivity to variation in environmental factors (which remain within optimal values for that species) is virtually absent. Thus, the specificity of the action of ecological factors simultaneously operating at optimal levels for a species disappears. In fact, the prevailing gradient of environmental conditions represents such a synthetic construct. The specific effect of a given factor can manifest itself when its values deviate from the optimum, falling into excess or deficit. Is the complete correspondence of species and their communities to the prevailing synthetic gradient of environmental conditions, therefore, a form of apparent insensitivity to the environment? The “pure” effects of ecological factors, when controlling for the covariation of other factors and the affiliation of communities to specific associations, were statistically significant. However, they explained only a small fraction of community inertia. The pure influence of factors has a pronounced fine-scale spatial component. In contrast, the action of “primary” factors in the context of synthetic factors operates at broad or meso-scale spatial levels. Thus, local mismatches between the ecological optima of plant species and the prevailing environmental conditions may account for the specific manifestation of the effects of primary ecological factors.

5. Conclusions

The core of phytoindication lies in assessing the response of biotic systems to environmental factors. Phytoindication is not a simplified method of measuring environmental properties that is less accurate than physical or chemical techniques; rather, it is a unique ecological tool for studying the impact of the environment on ecosystems. This is emphasised by the two main patterns of indicator values identified within the studied landscape system. The first pattern reflects the coordinated dynamics of natural environmental factors: an increase in soil moisture and nutrient availability correlates with a decrease in light, temperature, continentality, and soil acidity. The second pattern indicates variability caused by anthropogenic impact, expressed as a gradient from increased hemeroby to greater naturalness of plant communities. The established correlations between indicator values are not methodological artefacts but natural consequences of reducing the multidimensional space of factors to the response space of biotic systems, which integrate the effects of numerous factors. A specific response of plant communities to individual primary factors occurs only when these factors deviate significantly from the communities’ optimal values. This reaction manifests at finer temporal and spatial scales, where the effects of factors are no longer ‘masked’ by the integrative action of the biosystem. It is advisable to consider interconnection patterns in responses to environmental factors when comparing landscape systems, as absolute phytoindication scores are not directly comparable due to their interdependence.