The study site is an estuarine karst environment that can be characterized by wetland features. It is located inside Raša Bay (Figure 1) belonging to the north part of the Adriatic Sea, in a Štalije settlement, coordinates of 45°02′50′′N and 14°02′82′′E. The area is characterized by hard carbonate bedrock overlain by reddish clay-loam cambic type of soil. The karst topography and residue of variable thickness from the weathering of carbonate strata typify the terrain. The study area has a Mediterranean type of climate, with mild humid winters (annual precipitation of 900–1000 mm), hot dry summers, and prevailing NE winds (the so-called bora) (Durn et al., 1999; Peh et al., 2010). Since the local area is composed mainly of highly vulnerable karst, its environmental quality should be carefully monitored (Medunić et al., 2020a).

The sampling campaign was carried out on October 5, 2021. Sediment and soil samples were collected from sites (No. 1–8) within some 500–700 m from the RCSW as indicated in Figure 1. Samples No. 1-7 were collected from surface soil and sediment (0-5 cm). The position of samples No. 6–8 corresponds to the site of former lagoons (a naturally shallow part of the estuary) that received coal sludge (coal processing wastewater) some 60-90 years ago. There, huge quantities of wastewater were stored and gradually filled up with coal and mud particles. Today, it is partly covered by grass and shrub vegetation. Samples No. 2–5 were collected downstream from the lagoon along Raša River. Sample No. 1 was collected from the estuary sediment located halfway between the local road and the coastline, partly wetland. Three samples with different depths were collected at site No. 8 with a stainless soil core sampler that was partly filled with groundwater. The samples were designated as No. 8-10, 8-20 and 8-30, corresponding to the depths of 15–20, 20–30, and 30–40 cm. Following the return to the lab, the samples were divided into two portions. One subset was stored at −80 °C for microbial DNA extraction. The other subset was air-dried, and then carefully packed in sterile bags for further analyses.

The physiochemical properties of the collected samples and their PAH concentrations are shown in Table 1. The samples had generally alkaline pH (7.6-8.3). Depending on the collection sites, moisture ranged widely from 19% to 77%. Total organic carbon ranging from 98 to 268 mg/kg and C:N ratio between 4.8 and 26.0 were recorded in the samples. Elevated levels of TS (up to 9.3%) and OS (up to 4.2%) were detected in samples No. 6-8 (Figure 1), which represent the former coal sludge-receiving site. Interestingly, the pH, moisture, C:N ratio, TS and OS of sampling sites No. 1-5 were significantly different from No. 6-8 (ANOVA, P < 0.05).

Two types of PAHs including naphthalene and its derivatives [naphthalenes], and dibenzothiophene [DBT]) were quantified (Table 1). PAHs and heavy metals in an uncontaminated soil sample in the same area from our previous studies were served as references (Medunić et al., 2016). Results revealed that the samples contained extremely high concentrations of naphthalene (26.3-40.1 mg/kg) and DBT (7.4-8.9 μg/kg) from the coal sludge-receiving site compared to downstream locations (0 μg/kg) and control sample (0.1 μg/kg). The correlations of physicochemical properties with PAHs (Supplementary Table 1) were not significantly affected by the depth at site No. 8 except for TS and OS. In addition, for samples from all sites, correlation analysis revealed that C:N ratio and PAHs (Spearman ρ = 0.692, P = 0.027), TS and PAHs (Spearman ρ = 0.925, P = 0.000), OS and PAHs (Spearman ρ = 0.937, P = 0.000) were positively correlated (Supplementary Table 2). Although the concentrations of specific PAH congeners varied from site to site, the anomalously high levels of PAHs in No. 6-8 samples indicate that the site is still heavily polluted from the former coal sludge disposal.

Table 1 also shows that the concentrations of major HTEs in the samples are considerably higher compared with the references (Medunić et al., 2016). The concentrations of selenium, uranium and vanadium in samples No. 6-8 are significantly higher than that in samples No. 1-5, while lead and zinc showed the opposite (ANOVA, ρ < 0.05). Furthermore, positive correlations of Se-U-V (Spearman ρ > 0.69,P < 0.05) and Pb-Zn-Cu (Spearman ρ > 0.71,P < 0.05) were observed which might be worthy of further investigation.

A total of 610,122, 682,950, and 788,679 high-quality bacterial, archaeal, and fungal sequences were obtained and grouped into 8,231, 607 and 4,722 OTUs using a 97% sequence similarity cutoff. Then OTU tables were, respectively, rarified to the minimum depth of sequences of 30,709, 43,999, and 52,799 per sample.

Figure 2A shows alpha-diversity levels of bacteria, archaea, and fungi which were all higher in samples No. 1-5 than in samples No. 6-8, except for the Shannon index for the archaea. The community diversity and abundance of bacteria, archaea and fungi were not significantly affected by the depth, suggesting that bacterial, archaeal and fungal communities were not sensitive to the soil and sediment depth at site No. 8. These were confirmed by the correlation analysis (Supplementary Table 1). Furthermore, the differences in the alpha-diversity indices between the coal sludge-receiving site and the downstream locations were confirmed via the Wilcoxon rank sum test (Figure 2A). It appears that except for archaeal diversity, the microbial diversities and abundance were significantly higher in the downstream locations than in the contaminated site.

Non-metric multidimensional scaling (NMDS) analysis revealed that samples No. 1-5 and No. 6-8 formed distinct clusters (Figure 2B), with significant differences being found at OTU levels (analysis of similarity test [ANOSIM]). These differences among samples No. 1-5 and No. 6-8 were the highest for bacterial communities, followed by archaeal and fungal communities, indicating that the coal wastewater discharge impacted the bacterial communities the most.

Figure 3 shows the microbial community composition at the phylum level for a relative abundance greater than 0.1%. Bacterial sequences were primarily composed of the phyla Proteobacteria (27.92%), Chloroflexi (15.47%), Actinobacteriota (9.59%), Desulfobacterota (9.15%), Acidobacteriota (7.25%), and Bacteroidota (7.12%). The majority of the archaeal sequences belonged to the phyla Crenarchaeota (54.19%), Halobacterota (21.27%), and Thermoplasmatota (11.48%). The most abundant fungal phyla were Ascomycota (72.06%), Rozellomycota (7.45%) and Basidiomycota (6.98%).

Microbial community profiles at the genus level were demonstrated by a hierarchically clustered heatmap (Figures 4A–C). The top 30 genera and 10 samples were both hierarchically clustered based on the Bray-Curtis similarity index. Bray-Curtis dissimilarity analysis revealed that microbial communities of the 10 samples were clustered into two groups. The bacterial and archaeal clustering indicated that samples 1-5 (the downstream locations) were closely related while samples 6-8 (the coal sludge-receiving site) shared a close relationship. Conversely, samples No. 6, 8-30 cm, and 8-40 cm exhibited a distant relationship with samples No. 7 and 8-20 cm in fungal communities. This affinity relationship was also intuitively demonstrated by the NMDS analysis. As shown in Figures A–C, the most representative bacterial and archaeal genera in samples No. 6-8 were Blastococcus and Halogranum, respectively. In samples No. 1-5, Sulfurovum, Sulfurimonas, norank_f_Desulfurococcaceae, Methanosarcina, and Methanobacterium were the most abundant genera. Although most of the dominant genera of samples No. 1-5 and No. 6-8 were different, the top 30 genera were shared by all the samples. This may represent the ecological coherence of the water and sediment of Raša river.

Wilcoxon Rank Sum tests were conducted (double tail test, fdr multiple correction) to determine the difference in abundance of the classified bacterial, archaeal and fungal taxa between the sample groups of No. 1-5 and No. 6-8. The species with significant differences in the sum of the means in the top 15 genera are presented in Figures 4D–F. Interestingly, most of these bacteria and archaea were significantly enriched in samples No. 6-8. Specifically, Blatococcus, “Subgroup 10” from Thermoanaerobaculaceae, Nitrosococcus, Ferruginivarius, Nitrosopumilus, Methanomassiliicoccus, Halogranum, Halobacterales (order), Haladaptatus, Haloterrigena, and Natronomonas were enriched in samples No. 6-8 (all P < 0.05). The same members in samples 1-5 were all less than 0.1%. Conversely, most fungi were significantly enriched in samples No. 1-5, while only one genus showed an abundant advantage in samples No. 6-8, i.e., Neocamarosporium. In comparison, members of Sulfurovum, Anaerolineaceae (family), Rhodobacteraceae (family) and Bathyarchaeia (class) were enriched in samples No. 1-5. The members belonging to Sulfurovum and Bathyarchaeia (class) accounted for 3.54% and 60.6% of the population in samples No. 1-5 but less than 0.01% samples No. 6-8.

Canonical correspondence analysis (CCA) was used to determine the relationship between biochemical factors and bacterial, archaeal or fungal communities. Because the strong multicollinearity was detected among nine HTEs (variance inflation factor (VIF) value > 10), VIF analysis was performed prior to the CCA analysis (threshold value = 2). As shown in Supplementary Table 3, two HTEs’ environmental factors, Zn and U, were ruled out. As shown in Figure 5, for the eight environmental variables, acute angles were formed among variables of PAHs, U, C:N ratio (C/N), pH, and OS, indicating that these variables had a synergetic impact on the microbial community. On the contrary, sample moisture and Zn showed an opposite relationship in obtuse angles on the microbial composition with the other five variables (PAHs, U, C/N, pH, OS). Apparently, the microbial composition of samples No. 6-8 is positively correlated with the levels of PAHs, U, C/N, pH, and OS, and negatively associated with the sample moisture, thus implying that a higher level of PAHs, U, C/N, pH, and OS, and lower moisture might reduce the abundance and diversity of microbial communities in soils and sediments. All the edaphic variables explained 92.3, 93.4, and 88.1% of the bacterial, archaeal, and fungal variance, respectively. The key environmental factors impacting microbial communities in the samples were pH, OS, and PAHs (P < 0.05).

Provided that the microbial diversity and community composition are different between sample groups No. 1-5 and No. 6-8, the key functional microorganisms between the two groups were compared and presented in Figure 4G. Notably, 21 genes (including pcaG, pcaH, ligA, ligB, phdF), 23 genes (including adh, frmA, ADH5, nahF), and 16 genes (including dmsB, dmsC, ssuE, msuE, and cysJ) were affiliated with important enzyme categories (dehydrogenase, hydroxylase, aldolase, decarboxylase, hydrolase) involved in the degradation of PAHs (Supplementary Figure 1), naphthalene (Supplementary Figure 2), and sulfur-containing compounds (Supplementary Figure 3). These pathways were predicted to be harbored by the soil and sediment bacteria and archaea. Among these genes, the sulfur degradation-related genes generally showed the highest abundance followed by naphthalene and PAHs (Figure 4G). The abundance of functional genes associated with naphthalene and sulfur degradation in samples No. 6-8 was greater than that in samples No. 1-5. Notably, most of the genera with significant differences in the top 15 bacterial and archaeal genera have functional genes that are related to the degradation of naphthalene, PAHs, or sulfur-containing compounds. Specifically, the genus, Ferruginivarius can degrade all three substances while the other 5 genera (Natronomona, Blastococcus, Halogranum, Haladaptatus, Haloterrigena) could degrade two substances (Naphthalene/PAHs, sulfur-containing compounds). Both groups exhibited a large number of genes that are involved in sulfur metabolism. In addition, total enzymes and enzymes capable of PAHs degradation for fungi in the two sample groups have shown almost no difference (Table 2), suggesting that indigenous fungi contribute little to the degradation of PAHs in the long run.

Elevated concentrations of HTEs were documented in all collected samples (Table 1), thus indicating that the Raša Bay area in the vicinity of the RCSW is severely polluted with HTEs as a result of the legacy of coal sludge. Previous publications report that the local environment (e.g., groundwater, soil, wild plants, crops, and aquatic sediment) is polluted with HTEs due to former coal mining, processing, and combustion activities (Sando et al., 2015; Medunić et al., 2016, 2020a,2020b,2021). Interestingly, the concentrations of Mo, Cu, Zn, Cr, and Pb in samples No. 1-5, especially Zn and Pb, were higher than that in samples No. 6-8, suggesting migration occurred possibly from the erosion and upstream leaching of the coal sludge-receiving site (No. 6-8), while Se, U, and V are coalphilic elements (Yudovich and Ketris, 2006). However, Mo, Cu, Zn, Cr, and Pb are less coalphilic, and more soluble in water (Medunić et al., 2016). Therefore, leaching of Mo, Cu, Zn, Cr, and Pb from coal debris contained in the sludge-receiving site soil could have happened. A previous study (Medunić et al., 2016) also indicated positive correlation (P < 0.05) between Cu-Pb-Zn and Se-U-V clusters in soil polluted by Raša coal.

PAHs are components of fossil fuels such as coal. DBT and its derivatives are the main sulfur-containing aromatic compounds in fuels (Kertesz and Wietek, 2001). These aromatic compounds and their degradation products have been found in groundwater, seawater, sediment, soil and atmospheric samples related to old gas plants, industrial waste sites and wood treatment facilities (Warshawsky, 1992; He et al., 2009), which can persist in the environment (Xu et al., 2006; He et al., 2009). According to the classification proposed by Maliszewska-Kordybach (1996), the former coal sludge-receiving site (samples No. 6-8) is heavily polluted (>> 1,000 μg/kg, Table 1). On the contrary, PAHs were not detected in the downstream samples (samples No. 1-5). Notably, the PAH composition of coal is commonly dominated by an average of 4.5-ring PAHs (Dong et al., 2012), while those detected in the analyzed samples were dominated by a lower number of rings (2-rings). In addition, concentrations of PAHs were also lower compared with raw Raša coal (53,407.7 μg/kg) (Medunić et al., 2016). These results indicate that PAHs were not migration-prone and were naturally attenuated by the soil and sediment microbes (Schneider et al., 1996; Coates et al., 1997; Yagi et al., 2010). However, the concentration of naphthalene from the polluted coal sludge-receiving site was apparently higher than coal-tar waste (Yagi et al., 2010). This is likely due to slow desorption rates for PAH compounds in coal-contaminated floodplain soils (Achten et al., 2011).

Complex variation occurred in sample properties due to coal processing wastewater, such as increases in sample pH and C:N ratio. These changes could be partly explained by microorganisms, whose activities have proven essential for the functioning of these nutrient cycles (Falkowski et al., 2008; Dai et al., 2018; Jiao et al., 2018). Studies have shown that when pH was between 8 and 9, the activities of phenanthrene, pyrene and naphthalene-degrading microorganisms were highest (Song et al., 2012; Lian et al., 2015). Furthermore, a suitable C:N ratio is beneficial for the continuous growth and proliferation of PAH-degrading microorganisms in soil and sediment (Cai et al., 2021). However, these environmental changes can also lead to an alteration of microbial communities (Schneider et al., 1996; Chaput et al., 2015; Jiao et al., 2018; Liu et al., 2018), as shown by the canonical correspondence analysis (Figure 5). On the other hand, a high concentration of sulfur was detected in the polluted coal sludge receiving site, attributed to Raša coal which has exceptionally high values of sulfur (up to 9-11%). Noteworthy, previous study also confirmed that sulfur and PAHs were strongly correlated in soil polluted with Raša coal as a consequence of a local coal-fired power plant (Spearman, P < 0.01) (Medunić et al., 2016).

HTEs have been repeatedly illustrated to have significant impacts on microbial community structure and function across diverse habitat types (Nicomrat et al., 2006; Weber et al., 2008; Khan et al., 2010; Liu et al., 2018; Zhao et al., 2021). Although this study showed that the levels of HTEs were elevated in all samples, CCA analysis revealed that HTEs may not have played an important role in the change of microbial communities (Figure 5, all P > 0.05). On the contrary, the accumulation of PAHs has been the main cause for changing the structure of microbial communities and thus their functions. Since HTEs cannot be degraded like TOPs, they may not participate in the basic element cycle of cells (Nies, 1999). Furthermore, it has been reported that microbial cells have ubiquitously defensive mechanisms against HTEs in such environments (Nies, 2003). In this study, there are bacteria and archaea with intrinsic high-level resistance to HTEs that have been reported in other studies [e.g., Bacillus sp. (Amoozegar et al., 2005), Marinobacter sp., Halococcus sp. (Voica et al., 2016)].

This study shows that PAHs decreased the biodiversity of sample microbiomes, which is in accordance with the fact that high concentrations of PAHs negatively affect the total microbial populations (Zoppini et al., 2016; Picariello et al., 2020). PAHs have detrimental impacts on microbial metabolic processes, enzyme activities (Picariello et al., 2020), and cell functional integrity (Zoppini et al., 2016). For example, PAHs can inhibit oxidative phosphorylation processes by preventing the ATP-forming mechanism (Zoppini et al., 2016). This would explain why the PAH pollution of the Raša coal sludge-receiving site is one of the key factors affecting the microbial community at this study locality. Studies showed that fungal biomass and enzyme activities associated with PAH degradation (e.g., laccase, peroxidase, cytochrome P450 monooxygenase) increased in the early PAH contamination (Marco-Urrea et al., 2015; Mineki et al., 2015; Jové et al., 2016; Bellino et al., 2019; Picariello et al., 2020). However, a minuscule difference in enzyme activity related to PAH degradation was observed in this study. This suggests that fungi have few roles to play in the later stage of the degradation of PAHs, while bacterial communities were more responsive and adaptive to contamination than fungi. This could be ascribed to the high intrinsic growth rates and resilience in disturbances and perturbations of bacteria, rendering a more rapid response to environmental changes such as coal processing wastewater pollution in this study. This contrasting pattern could reflect that bacteria and archaea were more sensitive to PAH pollution from the coal sludge-receiving site (No. 6-8).

Recent research has revealed significant influences of coal waste-associated pollution on the composition and diversity of indigenous soil microbial communities (Ghiorse et al., 1995; Yagi et al., 2010; Roth et al., 2019; Hamidović et al., 2020). This study shows that Proteobacteria, Chloroflexi (bacteria), Crenarchaeota, Halobacterota (archaea), Ascomycota, and Basidiomycota (fungi) were the dominant phyla in the coal sludge-receiving site (affected by the RCSW). Similar phyla have previously been found to dominate the soil microbiota in coal-related wastes (Yagi et al., 2010).

Noteworthy, while the top 30 genera were shared by all samples, a very different distribution pattern of dominant genera was observed. Specifically, Blastococcus, “Subgroup 10” from Thermoanaerobaculaceae, Nitrosococcus, Ferruginivarius (bacteria), Nitrosopumilus, Methanomassiliicoccus, Halodaptatus, Halobacterales (order) (archaea), and Neocamarosporium (fungi) were significantly more abundant in the polluted coal sludge samples (P < 0.01). Interestingly, Blastococcus, Ferruginivarius, “Subgroup 10” from Thermoanaerobaculaceae, Methanomassiliicoccus, Halodaptatus, and Halobacterales (order) have shown the ability to catabolize complex aromatic compounds or simple hydrocarbons (Joo and Kim, 2005; Youssef et al., 2012; Dyksma et al., 2016; Zárate et al., 2021; Vidal-Verdú et al., 2022). Additionally, “Subgroup 10” from Thermoanaerobaculaceae (Basen et al., 2018), Nitrosococcus (Dyksma et al., 2016), Ferruginivarius (Wiese et al., 2020), Halodaptatus and Halobacterales (order) (Wang et al., 2022) are involved in the sulfur oxidation process. “Subgroup 10” from Thermoanaerobaculaceae has strong viability in polluted coal sludge soil and may have participated at the end of the sulfur cycle in the degradation of TOPs and benefit the growth of other microorganisms (Basen et al., 2018). Besides, Nitrosococcus can release nitrogen and nitrogen oxides (Watson, 1965), which explains why samples No. 6-8 have a higher C:N ratio. Subsequent functional inference of communities further revealed that Blastococcus, Ferruginivarius, “Subgroup 10” from Thermoanaerobaculaceae, Methanomassiliicoccus, Halodaptatus, and Halobacterales (order) may have relevant genes involved in the PAH and naphthalene degradation. “Subgroup 10” from Thermoanaerobaculaceae, Nitrosococcus, Ferruginivarius, Halodaptatus and Halobacterales (descending order) may have relevant genes involved in the sulfur-metabolism pathways. Overall, our results demonstrate that the native community in the study area contains microbes that can tolerate high concentrations of HTEs and degrade PAHs and sulfur-containing organic compounds. Notably, the microorganisms that can degrade refractory and sulfur-containing organics may have been enriched as they adapt to the surrounding environment. However, the restoration (biodegradation process) of the polluted ecosystem is slow as extremely high levels of PAHs and HTEs are still observed today, following some 60 years of natural attenuation processes. This study shows that polluting organic compounds mostly shape the microbial structures in the local estuarine environment.