In recent decades, engineered ecosystems characterized by considerable physical alterations of the local environments have attracted significant attention (Jones et al., 1994). Plantation ecosystems, including planted forests, are often viewed as a compensatory response to the continuous long-term decline of the total area covered by natural biotopes, mainly due to human interventions (Tilman and Lehman, 2001). However, typically limited plant composition represents their common drawback, leading to the reduction in the local soil microbiota diversity and stability (Chen et al., 2019), that in turn appears an essential driver of the flora diversification and vegetation productivity (Van Der Heijden et al., 2008). These complex symbiotic interactions largely govern stress resilience and adaptation capabilities of the entire ecosystem (Cardinale et al., 2012; Liang et al., 2016), especially under combined stress factors, imposing major challenges for the engineered ecosystems management. Therefore, a better understanding of the complex interplay between soil geochemistry, microbial biodiversity and flora vegetation dynamics that are inevitable components of the overall ecosystem stability and its adaptation capabilities is essential for overcoming the above challenges.

While the main focus of interest has been traditionally shifted towards anthropogenic effects, multiple studies of zoogenic and ornithogenic ecosystems largely focused on the unique and vulnerable habitats in the polar regions and on remote islands (Zmudczyńska-Skarbek et al., 2015; Zeglin et al., 2016; Rowe et al., 2017; Zmudczyńska-Skarbek and Balazy, 2017; Abakumov, 2018; Guo et al., 2018; Potapowicz et al., 2020; Abakumov et al., 2021; Nizamutdinov et al., 2021). Nevertheless, terrestrial ornithogenic ecosystems are also of considerable interest, since they often represent long-term natural testbeds of locally altered environments, and thus provide significant information on the adaptability of local flora and fauna, that can be studied without introducing any additional anthropogenic invasion.

Large bird colonies occupying relatively compact areas deposit excessive amounts of nutrients (Frederick and Powell, 1994; Ashworth et al., 2020), leading to considerable alterations in the soil geochemistry and microbiota, affecting the surrounding vegetation, and thus altogether considerably altering the local ecosystems (Sekercioglu, 2006; Whelan et al., 2015; Natusch et al., 2017; Lowney and Thomson, 2021; Grant et al., 2022; Hawke, 2022; Lowney and Thomson, 2022). For example, the impact of large birds like cormorants (Phalacrocorax carbo), Adélie penguins (Pygoscelis adeliae), Eurasian Crane (Grus grus) and Grey Heron (Ardea cinerea) has been intensively studied (Goc et al., 2005; Hobara et al., 2005; Kameda et al., 2006; Adamonytė et al., 2013; Klimaszyk et al., 2015; Klimaszyk and Rzymski, 2016; Guo et al., 2018; Matulevičiūtė et al., 2018; Veum et al., 2019; Al Shehhi and Muzaffar, 2021; Machač et al., 2022; Valkó et al., 2022). These birds form large colonies on coasts and forests leading to an intensive deposition of allochthonous substances in the local environment and consequent eutrophication, in turn altering the soil biogeochemistry, degrading the biodiversity and suppressing plants vegetation (Ishida, 1996; Anderson and Polis, 1999; Hobara et al., 2005; Kolb et al., 2012; Adamonytė et al., 2013). In turn, these colonies could be viewed as long-term natural testbeds often contributing to the local environment over decades, where the endpoint of this multi-year natural experiment could be observed and analyzed here and now.

Soil microbiomes have been proposed as sensitive indicators of global change and an integral part of biogeochemistry (Oyugi et al., 2006; Varin et al., 2012; Santamans et al., 2017). The structure of microbiome significantly depends on carbon and nitrogen (Ligeza and Smal, 2003; Harrow et al., 2006; Schaefer et al., 2014; Tytgat et al., 2016; Otero et al., 2018), moisture (Lavian et al., 2001), and phosphorus (Chong et al., 2009; Roesch et al., 2012; Kim et al., 2014). In turn, the forest floor quality and biogeochemistry deposition are among common consequences of the ornithogenic ecosystem alterations (Hobara et al., 2001; Ellis, 2005). Thus, the birds feaces, which are extremely rich in phosphorus (Osono, 2012; Zhu et al., 2014; Wurster et al., 2015; Domínguez et al., 2017; Telesford-Checkley et al., 2017; Otero et al., 2018), nitrogen (Legrand et al., 1998; Tomassen et al., 2005; Barrett et al., 2006; Aislabie et al., 2008) and ammonia (Mizutani and Wada, 1988; Zhu et al., 2011; Riddick et al., 2014; Crittenden et al., 2015; Croft et al., 2016), as well as organic matter (Huang et al., 2014; Chen et al., 2020), affect the soil and sediment microbiomes (Zhu et al., 2015; Shen et al., 2023). Moreover, they introduce significant concentrations of Mg, Ca, K and Zn (Ellis et al., 2006; Breuning-Madsen et al., 2010; García et al., 2011) into the affected soil, having significant impact on both microbial communities (Wang et al., 2015; Santamans et al., 2017; Minkina et al., 2022) and fungal diversity (Adamonytė et al., 2013; Kutorga et al., 2013) of the latter. These alterations initiate consequent considerable shifts in the soil macrofauna (Korobushkin and Saifutdinov, 2019) and formation of unique ornithophilic and ornithochorous vegetation (Bradbury et al., 2005; Azpiroz and Blake, 2016).

Conifer trees including Scots pines (Pinus sylvestris L.) are among the common nesting sites and thus are considerably affected by bird colonies (Kutorga et al., 2013; Żółkoś et al., 2013). While some investigations reported much better growth of pines at the sites with bird droppings (Tomassen et al., 2005), the long-term impact of bird feaces leads to the suppression of the trees and forest-specific plants growth, followed by their eventual replacement with ruderal and light-preferring plants due to both drastic changes in the soil biogeochemistry, as well as toxic properties of the bird feaces (Ishida, 1996; García et al., 2011). Finally, intensive vegetation of nitrophylic plants leads to the accumulation of dead biomass and diversity degradation (Jones et al., 1994; Mun, 1997; Moore, 2006).

Despite a large number of studies on the impact of bird colonies on soils and surrounding vegetation, many of them have limited focus either on soil geochemistry or on its microbiota or vegetation dynamics. In this work, we aimed at a better understanding of the complex interplay of various geographical, geochemical, microbiological, and dendrological characteristics, as well as their manifestation in the vegetation indices obtained by remote sensing observations in a Bayesian network statistical model representing the overall impact on the local environment.

We consider a regularly planted pine forest area as a prominent testbed, where the emergence and the following 15-year long expansion of a Grey Heron colony could be viewed as a nearly perfect in situ natural experiment, with surrounding unaffected forest of initially the same age, composition and diversity representing a relevant control. We analyze how long-term deposition of nutrients from the bird colony altered the local soil biogeochemistry both under the current nesting sites and in the surrounding area, leading to major alterations in the local microbiota. In turn, these conditions affected the trees and other surrounding flora growth, that has been particularly reflected in their increased vulnerability to combined biogeochemical and climate stresses, as well as reduced adaptational capabilities indicated by their limited recovery from the major 2010 drought. Finally, we also propose a statistical model reflecting the complex interplay between various factors from initial biogeochemical contributions by the bird colony to their effect on the vegetation dynamics and its reflection in the multispectral vegetation indices, and show how this model could be potentially employed in a couple of prediction scenarios.

Increasing anthropogenic influence on the environment drastically changes the ecological profile of areas with high population density and/or intensive agricultural activity, that in turn requires real-time evaluation and prediction of the ecosystem status for a timely risk assessment and informed environmental management. While the engineered ecosystems partially compensate the decline of the natural biotopes, on the other hand, they are characterized by considerable alterations in the local environments and soil biogeochemistry (Jones et al., 1994; Tilman and Lehman, 2001).

In this research, we analyzed complex interrelations between the soil geochemistry, microbial biodiversity and flora vegetation dynamics in a planted Scots pine forest hosting a large colony of Grey Heron, see Figure 1 for location details. Being nested from approximately 2006, the colony consists of approximately 280 bird nests (as of summer 2022) located on the tops of mature Scots pines trees being surrounded by a relatively unaffected forest area. Given the initially uniform characteristics of the regularly planted forest area, the surrounding forest of the same age, composition and diversity, but barely affected by the bird colony, represents a relevant control. Thus, the investigated Grey Heron colony appears a prominent natural in situ experiment unfolding the changes in an engineered ecosystem in response to the excessive deposition of both biogenic elements (N, P) and toxic organics in the soil, as well as considerably reduced pH values compared to the presumably unaffected surrounding forest area. The above conditions are representative for both natural areas occupied with birds (Goc et al., 2005; Hobara et al., 2005; Kameda et al., 2006; Adamonytė et al., 2013; Klimaszyk et al., 2015; Klimaszyk and Rzymski, 2016; Guo et al., 2018; Matulevičiūtė et al., 2018; Veum et al., 2019; Al Shehhi and Muzaffar, 2021; Machač et al., 2022; Valkó et al., 2022) as well as agricultural environments (Bradbury et al., 2005; Minkina et al., 2022). Accordingly, considering the emergence of the bird colony around 2006 and its further expansion as an experimental input, one can currently observe an engineered ecosystem with accumulated alterations in soil biogeochemistry and consequent impact on the trees and other surrounding vegetation as an endpoint of a 15-year long natural experiment.

In order to select the keynote factors that govern the ecosystem changes, and thus could be used as indicators for predicting significant alterations in the environment, a cross-correlation matrix between them was obtained (Figure 6), that in turn allowed to design a non-parametric Bayesian network model (Figure 7) revealing the complex interplay between various ecosystem components, where the nodes are associated with particular geochemical, microbiological, dendrological and vegetation indicators, while the edges are characterized by their interrelations quantified by respective rank correlations R. In order to simplify the model by eliminating weak and often secondary effects, in its graphical representation we show only those correlations R exceeding an empirically chosen threshold 0.6.

Degradation of biodiversity and suppression of plants vegetation in areas occupied with large bird colonies has been reported by many authors (Hobara et al., 2001; 2005; Adamonytė et al., 2013). In the investigated Grey Heron colony, apparently because of high toxicity of the bird feces (Ishida, 1996; García et al., 2011) and extreme pH decrease to 4.5 (Figure 2), the death of trees under the nests was observed in the initial location of the colony with the maximum exposition time to excessive nutrients accumulation, followed by the formation of a consecutively widening circular area with the expansion of the bird colony on the surrounding trees (Figure 1C). This fact is also reflected in the positive correlations between the trees growth dynamics (TRW index) and pH (R = 0.65, Figure 6), indicating negative impact of the soil acidity on the pine trees growth. Notably, there are also positive correlations between the trees growth dynamics and microbial diversity represented by the Simpson’s index (R = 0.29), suggesting that the repression of trees growth is associated with the reduction of the soil microbiota, in agreement with recent literature data. As well, the overall vegetation characterized by the MTVI2 index is positively correlated with pH, while demonstrating a weaker association than for the pine trees growth (R = 0.31 vs. 0.57), since other types of vegetation appear to be less sensitive to soil acidity.

As well, because of the influence of the bird feces, feathers, food remains, etc., dramatic changes in the soil biogeochemistry have been observed. In particular, significant increase of biogenic elements like nitrogen (∼10 fold) and phosphorous (∼4 fold) in the soil could be observed, while the organic carbon (C_org) was even reduced compared to the inner and the outer areas presumably reflecting the repression of the vegetation (see Figure 2). Interestingly, while maximum concentrations of nitrogen have been observed in the area under the trees currently nested by the birds, phosphorous exhibited its maximum concentrations in the inner circle corresponding to the initial location of the colony, that could be attributed to its presumable sequestration in non-soluble organic matters, an effect that has been observed in earlier studies (Brenner et al., 2006; Turner et al., 2006).

Suppression of the vegetation, directly reflected in the measured vegetation indices MTVI2 and BAI obtained from multispectral remote sensing data, see Supplementary Figure S9 and Figure 4, is also in agreement with the tree-ring data indicating that the trees growth rates were considerably reduced in the inner circle and in the nesting areas. As one can see from Figure 5, there are pronounced discrepancies between the TRW dynamics in the affected and unaffected areas. Remarkably, in the year 2008 trees in the area affected by the bird colony exhibited higher growth rates than in the control areas, presumably due to the fertilizing effects of the bird manure containing excessive concentrations of nitrogen (N) and phosphorus (P). However, despite of the obvious discrepancy, the growth rate altogether reduced in both areas continuously from 2008 until 2010 that coincides with the gradual enhancement of summer temperatures reaching its apogee in a major heatwave in 2010 accompanied by a flash drought indicated by a gradual PDSI decline.

To further quantify the combined biogeochemical and climate stress resilience of the trees, we next consider the 10-year time span between 2010 and 2019 corresponding to the recovery from the prolonged drought prior and during the extreme 2010 heatwave. We find that about one-half of all TRW records in the control area indicated significantly positive trend, while only three records in the impacted area followed a significant recovery pattern, see Supplementary Figure S11B. Remarkably, five TRW records in the affected area demonstrated negative x estimates, including one statistically significant negative trend, indicating severely reduced ability to recover from the drought stress under concomitant impact from the biogeochemical alterations induced by the bird colony.

In addition, statistically significant differences between the inner circle (inner), the nesting area (nests), and the surrounding unaffected forest (outer) could be observed not only in the relative trend x, but also in the average TRW data over the same 15-year period, as well as over shorter time spans, in particular, since the point of divergence (2007–2019), as well as over a more recent 5-year time span (2015–2019), see Supplementary Figure S12. However, high correlations between all above metrics reveal their strong interconnection, and thus only the relative trend that exhibited the most pronounced and significant discrepancies between studied sites has been considered for further analysis (see Supplementary Figure S12).

Furthermore, despite of the excessive deposition of nitrogen and available phosphorous with bird feces (Osono, 2012; Wurster et al., 2015; Crittenden et al., 2015; Domínguez et al., 2017; Telesford-Checkley et al., 2017; Otero et al., 2018), reduced microbial diversity could be observed in the soils of the colony area (Supplementary Figure S3) (R = 0.45–0.46, see Figure 6). These experimental observations are reflected in the Bayesian network model reconstruction (Figure 7). The model reveals explicitly that the colony has a great impact on pH, nitrogen and phosphorous of the soil, which in turn affects all other indicators like trees growth rate (TRW), microbial diversity (Simpson’s index) and some given bacteria.

Besides of the extreme deposition of nitrogen and phosphorous, the bird colony areas are characterized with increased concentrations of K, Na, Mg, Ca, Fe and Zn (Ellis et al., 2006; Breuning-Madsen et al., 2010; García et al., 2011) in soils, that in turn also affects soil microbiota (Wang et al., 2015; Santamans et al., 2017; Minkina et al., 2022). In our investigation, the ∼1.2 fold increase in the concentration was observed only for K in the nested area (R = −0.44, although statistically insignificant) and Zn (R = −0.61) in the central area of the colony (Supplementary Figure S2). Therefore, these elements are not linked with the nests position in the network (Supplementary Figure S10), although they affect the microbiota structure. Of note, the concentration of Zn is linked with trees growth and vegetation indices BAI and MTVI2. Additionally, while the quantities of Li, Co. and Mn were significantly reduced in the inner area compared to the nests and outer areas, only Mn was found to significantly correlate with microbial diversity. Remarkably, the concentrations of Li, Co., Fe and Ca are positively correlated with the microbial diversity (R = 0.32…0.48), while concentations of Li and Co. decrease in the inner area of the colony (R = 0.33 and 0.52, respectively). Finally, geochemical anomalies are reflected in the BAI index that is positively associated with the concentrations of N, P, Zn and Fe (R = 0.27…0.62), while negatively correlated with pH (R = −0.43, although some associations appear relatively weak and thus are not shown in the network).

The relevance of Arsenic (As) and Samarium (Sm) in the nests area is questionable and their increase could be a secondary indicator of the accumulation of an unlocalized contamination. Arsenic has been reported earlier to accumulate in bird tissues and feces from the environment (Sánchez-Virosta et al., 2018; Eeva et al., 2020). Although Samarium is widely investigated for the development of both antimicrobial and antitumor agents (Morais et al., 2014; Kratochwil et al., 2021; Zahmatkesh et al., 2022), no data regarding its accumulation in bird tissues are available.

The above alterations in soil geochemistry, especially the dramatic contrast of pH (Green et al., 2012) accompanied by a major excess in nitrogen and phosphorous (Teixeira et al., 2013; Gao et al., 2021), apparently represented the main factors leading to the changes in soil microbiota structure, compared to presumably unaffected areas (Figure 7). The increased abundance of Proteobacteria and Firmicutes, as well as decrease in Acidobacteriota and Actinobacteriota fit with similar changes in soils occupied by bird colonies and increased nitrogen reported earlier (Teixeira et al., 2013; Gao et al., 2021), represented by terminal nodes in the reconstructed network (Figure 7). While the linkage of Verrucomicrobiota, Gemmatimonadota and Myxococcota with birds has not been described yet, the decrease of this normal soil flora in the colony area (Supplementary Figure S2) is significant (R = 0.34…0.59, Figure 6), and thus could reflect negative changes in geochemistry and structure of the soil, and could act as relevant biomarkers of the changes in the soil ecosystem (Supplementary Figure S8). Nevertheless, only Gemmatimonadota is included in the final network based on the desired correlation threshold.

Finally, four genera Chujaibacter, Rhodanobacter, Bacillus and uncultured bacteria of Gaiellales were the most abundant ones in soils under the nests (Figure 3). Among them, only Chujaibacter, Rhodanobacter, and Bacillus could serve as biomarkers, since they are linked with many other nodes in the network (Figure 7). Relatively high abundance of Bacillus in the soil could be a consequence of their introduction with bird feather (Sotnychuk et al., 2020), that is reflected in correlation with nests location (R = −0.56). By contrast, significant increase in the abundance of Chujaibacter, the bacterium reported to be efficient in sulfur transformation and carbohydrates degradation (Zhang et al., 2022), and Rhodanobacter, reported earlier as a resident of highly contaminated areas with low pH (Green et al., 2012) is observed in the inner area of colony (R = −0.63…–0.68, Figure 6). In turn, Gaiellales were shown to be dependent of organic and fatty acids content (Sun et al., 2022).

While detailed investigation of the interconnections among various factors governing the local environment is essential for the understanding of their complex interplay represented by the graph in Figure 7, it appears timely and costly. A more practical application scenario would be to employ the revealed association patterns for the risk assessment and prediction making in local environments based on selected measurements only. From the environmental management point of view, the typical question would be, how a given natural (such as an emergence of a bird colony) and/or anthropogenic (such as introduction of a bird farm or another agricultural activity) invasions could impact similar ecosystems? Answering these questions would require simulation of different scenarios leading to the projected concentrations of nutrients and other chemicals introduced by the invasion and estimating potential impacts on the soil microbial diversity and vegetation dynamics for these scenarios.

Figure 8 shows one test example containing a reduced cross-correlation matrix (panel A) and a corresponding NPBN graph (panel B) indicating the interrelation between the concentrations of the keynote chemicals in the soil, the diversity of the soil microbial community and the local tree growth dynamics. To test the potential predictability, we employed Bayesian inference with a bootstrap-based validation scheme according to (Koot et al., 2023) with the upper layer of the graph (pH and keynote chemicals) in Figure 8B representing the input variables, while predicting the variables in the lower layer of the same graph. Figure 8C shows the results of the regression analysis between the observational and the predicted soil microbial diversity and tree ring width trend metrics, both characterized by correlation coefficients R > 0.8. We believe that the above result indicates the capability of the putative evaluation and prediction of future trends in ecosystem adaptation and environmental alterations based on the proposed approach.

To summarize, this study investigated the complex impact of a major Grey Heron colony on a planted Scots pine forest viewing it as a natural in situ experiment on an engineered ecosystem with limited plant and microbial diversity, leading to its increased vulnerability to combined biogeochemical and climate stress. We have shown explicitly that the long-term deposition of nutrients from the bird colony starting from its initial location at the time of emergence around 2006 and further expansion over the surrounding area in the next 15 years considerably altered the local soil biogeochemistry. Among the most pronounced effects, pH reduced to 4.5, as well as nearly 10-fold higher concentrations of phosphorous and 2-fold higher concentrations of nitrogen in their maximum accumulation areas compared to the presumably unaffected surrounding forest area could be observed, while C_org did not change significantly. Moreover, contrasts in the soil microbial diversity with considerable discrepancies in the relative abundance of several bacterial genera (Chujaibacter, Rhodanobacter, Bacillus, uncultured bacteria of Gaiellales and Acidobacteriales) has been observed both in the inner area of the colony and under the nested trees, in turn altogether affecting the trees growth and other surrounding flora vegetation. The latter is reflected in MTVI2 and BAI vegetation indices obtained by remote sensing, and apparently unravels the increased vulnerability of the forest ecosystem to combined biogeochemical and climate stress, as indicated by the limited recovery of the affected trees from the major 2010 drought stress, in marked contrast to the surrounding forest in terms of the annual growth rate (p = 3∙10⁻⁵, Mann-Whitney U-test). Based on the results obtained, we proposed a combined non-parametric Bayesian network model reflecting the complex interplay between geographical, geochemical, microbiological and dendrological characteristics, as well as their manifestation in the vegetation indices obtained by remote sensing observations. Finally, based on the revealed interrelations, using the Bayesian inference approach with a bootstrap based validation scheme, we have shown that the NPBN is capable of predicting both Simpson’s diversity index and tree ring width trend from just a few input variables representing the concentrations of the keynote chemicals in the soil with correlation coefficient R between observations and model based predictions exceeding 0.8. The above result indicates the capability of the putative evaluation and prediction of environmental trends in the ecosystem status and could further facilitate an informed adjustment of the planted forest management by an early planning of possible intervention scenarios.