Metagenomic analysis of the microbial community structure in protected wetlands in the Maritza River Basin

Abstract

Microbial communities drive the biogeochemical cycles in wetlands and provide a number of ecosystem services. They underpin soil function, and are easily impacted by anthropogenic pressure. This study examined the bacterial microbiome in the natural wetland of Zlato Pole and the protected, periodically flooded rice paddies in the Maritsa River Basin using a metagenomic approach, based on high-throughput sequencing. Alpha-diversity analysis showed a significant variation between the three study sites for Chao1 and ACE (abundance based coverage estimator) richness estimators. A positive correlation was established with pH, with highest values detected for the rice paddies and the lowest, for the Zlato Pole sediments. The obtained sequences were assigned into 37 known bacterial phyla with over 97% bacterial sequences classified within only nine of them. The bacterial communities in rice paddies sediments were more evenly distributed, whereas the Zlato Pole sediment was the most biased. The consortiums in the rice paddies were dominated by Proteobacteria, followed by Actinobacteria and Acidobacteria. The bacterial assemblages in those sites could not be distinguished by analysis of similarity. The Zlato Pole sediment community held an isolated position, where Acidobacteria was replaced by Firmicutes and Bacteroidetes and the microbiome showed an extremely high abundance of Gammaproteobacteria and Bacilli. The dominance of Gammaproteobacteria and the presence of Deinococcus-Thermus phylum, along lower nutrient concentration and the absence of correlation with the environmental parameters, is evidence of constant anthropogenic pressure around the area. The obtained results could be applied in the plans for sustainable management of the protected wetlands.

Keywords:

• Metagenomics
• 16S rRNA
• wetland soil
• community structure
• microbial ecology

Introduction

Wetlands are unique ecosystems combining terrestrial and aquatic habitats. As such, they exhibit some of the characteristics of each system, covering 7.7% of the existing land area [1]. They are estimated to be one of the most important terrestrial ecosystems providing a vast diversity of services, such as removal of pathogens [2], environmental pollution and flood control, nutrient cycling [3], terrestrial runoff buffers, etc. [4]. Natural wetlands globally are extremely sensitive threatened resources due to the constant urbanization, agricultural pressure and climate change [5]. Bulgaria is one of European Union member states with the lowest share of land area covered by wetlands (0.8%) [1]. The Maritza River Basin is a major inland area including different aquatic environments, such as Zlato Pole natural wetland, situated in the old riverbed, and periodically flooded rice paddies. Both wetlands are included in the national action plan for the conservation of wetlands of high significance in Bulgaria as potential Ramsar sites [6]. They play a crucial role as wildlife habitats and fall under the protection of the Birds Directive 2009/147/EC [7]. Their nature and situation in densely populated areas render them extremely sensitive to water and forestry management activities, extraction of sand and gravel, poaching and pollution [6]. The Zlato Pole protected area is highly vulnerable, as most of its natural ecosystems are or have been subjected to long-standing and strong anthropogenic press. On the other hand, paddy fields make up the largest anthropogenic wetlands on earth [8]. The differences in the aquatic and exploitation regime in relatively closely situated wetlands make them suitable for research on the effects of anthropogenic pressure on wetland health, which was the topic of the present study.

Wetland performance is often evaluated by direct measurement of vegetation or macroinvertebrates, but they do not provide an accurate assessment [9]. In recent years studies have been focused on the structure and diversity of sediment microbial communities as biological indicators [4,10], as they have a significant relationship with trophic status [11]. The sediment communities are more suitable for ecological assessment due to the much higher taxa diversity and abundance compared to the corresponding water body [12]. This is a result of the presence of complex matrix providing a surface for the microbial growth [13]. The microbial consortium (microbiome) in the soil at a given location includes diverse species based on ecological selection principles [14]. Two major groups of factors regulate the microbial community composition and diversity: the “bottom-up controllers”, including environmental factors such as nitrogen forms, organic matter, phosphorus, etc., and the biological factors (viruses, Protista grazing), known as “top-down controllers” [15]. Therefore, the sediment community richness and abundance are a reflection of the total sum of factors affecting the given ecosystem, and understanding the differences in the taxonomic composition in different types of wetlands is crucial [11]. The sediment communities have a vital role in the biogeochemical processes, energy transfer, and climate regulation. They can shape the physical and chemical characteristics of a wetland soil [13,16], having both direct and indirect effects on the health of the ecosystems [14,17].

This study aimed to provide a comparative analysis of the bacterial communities in Zlato Pole wetland and in seasonally flooded paddy fields, using a metagenomic approach, based on high-throughput sequencing, and to examine the correlation of the microbial communities’ composition with the soil properties. The knowledge of the soil microbiome could be applied in the plans for restoration, protection and sustainable management of the protected wetlands.

Materials and methods

Experimental locations

Wetland sites are located along the Bulgarian part of Maritsa River and include Zlato Pole wetland (ZP) and Tsalapitsa rice paddies (“Orizishta Tsalapitsa”), both protected under the Birds Directive (79/409/EEC). Protected zone Zlato Pole was used as a reference site with shores covered primarily with Typha sp. and Salix sp. It is the largest natural wetland (BG0002103) along the Bulgarian part of the Maritsa River. The site includes the old riverbed and several basins located several kilometers away from the town of Dimitrovgrad, south of the Zlato pole village. Two isolated, periodically flooded, rice paddies were selected in the protected zone “Orizishta Tsalapitsa” (BG0002086). It is a complex of rice paddies situated between the village of Tsalapitsa and the city of Plovdiv. The sampling locations are coded as follows: ZP1, Zlato Pole sediments (42°2.207', 25°42.938'); ZP2, soil from the non-flooded part of Zlato Pole wetland (42°2.241', 25°42.944'); C1, Tsalapitsa rice paddy sediments(42°13.600', 24° 33.804'; C2, soil from the non-flooded part of Tsalapitsa rice paddy (42°13.396', 24° 33.202'); P1, Plovdiv rice paddy sediments (42°10.307', 24°40.502'); P2, soil from the non-flooded part of Plovdiv rice paddy (42°10.436', 24°40.403').

Soil sampling

Soil samples from the topsoil 0–10 cm were collected in triplicates in July 2018 during the rice maturity stage, from flooded and non-flooded rice paddies and from the sediments and the non-flooded area at Zlato Pole. Each sample was mixed from nine random soil cores [18]. All samples were homogenized and sieved through a 2-mm mesh and partitioned into two subsamples: one was partially air-dried for chemical analysis, and one was placed in sterile 50-mL containers and stored at 4 °C in the dark for no longer than 24 h for microbiological analysis.

Physicochemical analysis

Soil moisture was determined by measuring the weight loss on drying the sample at 105 °C for 24 h. Soil pH was measured with a pH meter (WTW/SET) in a slurry of soil and water 1:2.5. Nitrate (NO₃-N) and ammonium (NH₄-N) nitrogen content was analyzed after soil extraction with 0.1 N KCl according to Motsara and Roy [19]. Total nitrogen content was determined by the Kjeldahl method [20]. Organic matter (OM) content was measured by calculating the loss of weight on ignition at 600 °C. Available phosphorus (AP) content was determined spectrophotometrically by Olsen’s method [19].

DNA extraction

Total genomic DNA was directly extracted from 0.25 g soil sample with DNeasy PowerSoil Kit (Qiagen) according to the producer’s protocol. The quality and quantity of the randomly-sheared, high molecular weight metagenomic DNA were evaluated by BioSpectrometer (Eppendorf) with µCuvette G 1.0 (Eppendorf) and agarose gel electrophoresis.

16S rRNA amplicon sequencing and data analysis

The V3–V5 hypervariable region of the 16S rRNA gene was amplified and sequenced using MiSeq Illumina platform with the 2 × 300bp paired-end read at Eurofins MWG Operon (Ebersberg, Germany).

The quality of the raw data was checked using FastQC. To preserve only high-quality sequences, all reads in the FASTQ files with sequencing errors or with ambiguous bases (“N”) were removed. After preprocessing, sequences were trimmed using Trimmomatic to 255 bp length to remove low-quality bases from the 3′ end [21]. Chimeras were detected and removed using USEARCH [22] as it is implemented in the QIIME-pipeline [23]. Operational taxonomic units (OTUs) were picked using the QIIME script (pick_open_reference.py); sequences were clustered with UCLUST and taxonomies assigned based on the SILVA (v132, https://www.arb-silva.de. . . . . . . . ) database at 97% identity cut-off value [24]. The OTU counts were normalized by subsampling to the lowest number of OTUs found in the sample. The community alpha diversity indexes were calculated using QIIME v1.9.1 [23].

Interactive visualization charts of the taxonomic diversity of the samples were generated using KRONA tool [25].

Data availability

Raw sequences for each sample were deposited in GenBank SRA archive and are available with the following accession numbers: BioProject accession number PRJNA563557; Samples accession numbers: SAMN12676176, SAMN12676177, SAMN12676178, SAMN12676179, SAMN12676180, SAMN12676181 (http://www.ncbi.nlm.nih.gov/bioproject/563557. . . . . .).

Statistical analysis

Statistical analyses were done using Primer 6 (Primer-E, Ltd), Statistica v.10 (StatSoft) and Microsoft Excel with XLSTAT software package (Addinsoft). Prior to the analysis, physicochemical parameters were standardized by subtracting the distribution mean from each value and then dividing the values to the standard deviation. Principle Component Analysis (PCA) of normalized soil chemistry variables was chosen to show the environmental parameter variations between samples. For the beta-diversity analyses, singleton and doubleton sequences were removed and the individual taxa abundance was transformed by the square root to downweight high abundance sequences [26]. Permutational multivariate analysis of variance (perMANOVA) and analysis of similarities (ANOSIM) were used to test the differences between the two types of wetlands as well as between the flooded and non-flooded soils in the studied areas and to assess the within-group and between-group means [27]. Non-metric multidimensional scaling (NMDS) plots were generated to visualize the dissimilarity in the microbial communities between the different sites. The NMDS was based on a Bray-Curtis similarity matrix [28,29]. Redundancy analysis was used to establish the relationships between microbial communities’ composition, soil type and predictor variables (environmental variables). Correlation analysis (Pearson) was used to analyze the correlation between genera and physicochemical parameters. The significance of all statistical tests was assessed at α = 0.05. Heatmaps and cluster dendrograms based on Euclidean dissimilarity metrics were constructed based on the relative taxa abundance [30].

Results and discussion

Physicochemical parameters

The physicochemical characteristics of the sampled soils are shown in Table 1. Principal component analysis was performed to compare the physical and chemical variations of flooded and non-flooded soils in the studied sites. The patterns are shown in Figure 1. The first two principal components explained 79.54% of the total environmental data variation. All samples, with the exception of Zlato Pole sediments, were situated along the PC1 axis. The ANOSIM categorized the flooded and non-flooded samples into two distinct groups, which were confirmed by the high R-value (R = 0.926; p < 0.01), with no significant differences between the natural wetland and the paddy fields (p = 0.133).

Figure 1. PCA soil biplot based on the environmental parameters sampled from the flooded and non-flooded areas of the studied wetlands. ZP1, Zlato Pole wetland sediments; ZP2, soil from the non-flooded part of Zlato Pole wetland; C1, Tsalapitsa rice paddy sediments; C2, soil from the non-flooded part of Tsalapitsa rice paddy; P1, Plovdiv rice paddy sediments; P2, soil from the non-flooded part of Plovdiv rice paddy.

Table 1. Values of the environmental parameters from the sampling sites in the studied wetlands.

Station	Water (%)	pH (H2O)	NH4-N (μg. g^-1)	NO3-N (μg. g^-1)	TN (%)	PO4-P (μg. g^-1)	ОМ (%)
ZP1	36.26	7.09	18.8	7.1	1.36	62.5	2.43
ZP2	17.81	7.28	76.2	12.4	11.04	74.9	48.35
P1	44.22	7.06	53.6	10.1	5.24	192.4	10.18
P2	8.40	7.25	56.9	37.2	10.12	89.2	52.38
C1	40.47	5.77	50.0	8.0	3.76	186.1	9.08
C2	8.3	6.85	58.7	26.3	6.45	78.0	76.27

The soils from Zlato Pole were sandy in structure, whereas the samples taken from the paddy fields varied from clay (flooded soils) to silty-clay (non-flooded samples). Non-flooded soils were characterized by neutral pH values, whereas in the flooded samples, pH varied between 5.77 (C1) and 7.09 (ZP1). The non-flooded soils were significantly more abundant in OM, TN, NH₄-N and NO₃-N, whereas AP was higher in the flooded paddy fields. In contrast, the flooded area of Zlato pole wetland tended to have low levels of OM, TN, and AP. In general the studied wetland sediments were more typical of natural than of constructed wetlands [31].

Alpha diversity of the bacterial communities

The analysis based on MiSeq Illumina sequencing of the V3–V5 region of the 16S rRNA gene generated 238 239 reads form the flooded area samples and 202 901 reads from the adjacent non-flooded area samples, which were reduced with 3% after removing low-quality sequences. After filtering, and OTU picking process, 181 328 sequences were assigned to taxa from the flooded area samples and 158 260 from the non-flooded samples. The total number of obtained high-quality OTUs per sample was used to calculate the α-diversity indices and the rarefaction curves (Table 2). Alpha diversity analyses were further used to determine whether the observed trends were driven by the studied environmental variables. There were evident differences between the reference wetland (ZP) and the paddy fields in terms of Chao1, ACE and Shannon diversity estimators. The rarefaction curves for all indices showed that all samples reached a plateau (Supplemental Figure S1), suggesting that all the diversity was retrieved from the sites. The Tsalapitsa rice paddies soil samples (C) had the steepest rarefaction curves with the highest taxonomic richness (Supplemental Figure S1). The values of the α-diversity indices supported these results: the highest values of ACE and Chao1 indices were in the C2 soil sample and the lowest ones, in the Zlato Pole sediments. There was significant variation between the three sites, both for Chao1 (F = 13.7774; p = 0.0018) and ACE (F = 15.6812; p = 0.0012) richness estimators, but there were no significant differences between the flooded and non-flooded samples (Chao1: F = 1.1766; p = 0.3035; ACE: F = 1.1705; p = 0.3047). Further analysis based on the pairwise comparisons confirmed that the C site had significantly higher Chao1 (p ≤ 0.05) and ACE (p ≤ 0.05) estimates than the P and ZP regions. The Shannon diversity index also varied between wetlands (F = 5.0921; p = 0.0332), with the lowest values established in the ZP1 sample (p ≤ 0.05). The results for Chao1, ACE and Shannon diversity estimators in the rice paddy fields exceed the previously reported values for floodwaters in vegetative and reproductive rice stages [32]. These findings are in agreement with Barber et al. [33], who reported lower diversity in native soils compared to agriculture fields.

Table 2. Number of bacterial taxonomic units recovered from the studied wetlands and diversity indices of 16S rRNA gene relative to each sample.

Station	Observed OTUs	Division	Class	Family	Genus	ACE	Chao1	Shannon
C1	6682	31	112	423	856	8300.733	7855.868	11.083
C2	7166	26	101	397	745	8993.878	8510.868	11.297
P1	6039	29	107	371	689	7531.278	7201.193	10.629
P2	5605	27	94	354	687	7313.096	6925.022	10.688
ZP1	4474	22	75	256	501	5813.220	5483.552	7.937
ZP2	5945	26	93	302	596	7249.514	6922.677	10.580

The ACE, Chao1 and Shannon diversity indices showed correlations with some of the analyzed physicochemical parameters (Supplemental Table S1). There was a strong positive correlation of pH with Chao1 and Shannon (p ≤ 0.01), and a weakly positive correlation with the ACE estimator (p ≤ 0.05). All three indices were weakly negatively correlated with the concentration of AP. There were no significant relationships between α-diversity and organic matter or nitrogen form concentration. In this regard, the increase in α-diversity may be a result of rice cultivation practices, such as fertilization. The provision of otherwise scarce nutrients, like nitrogen and phosphorus, may lead to a fertilization-induced increase in diversity [34].

Beta diversity among the studied areas

A comparison of the community taxonomic patterns based on multivariate analysis showed clear differentiation between the sites. This was also confirmed by NMDS ordinations, which showed separation of the bacterial communities based on soil type and area location (Figure 2). The permutational MANOVA design test and ANOSIM R values, based on the soil type, wetland type and location, substantiated this result and showed significant differences in the taxa composition between the natural wetland and rice paddies (R = 0.771; p ≤ 0.001), and between the flooded and non-flooded soil samples in the rice paddies (R = 1; p ≤ 0.05). After 10 permutations of random subsampling at each wetland site, the ZP community was consistently found to be significantly distinct from the rice paddies. In contrast, the bacterial assemblages present in the two rice paddies could not be distinguished (ANOSIM, R = 0.065; p = 0.124). The level of dissimilarity of bacterial communities observed in selected areas indicates significant correlations to the biogeochemical parameters in the natural wetland and the paddy fields. Hierarchical clustering based on the taxonomic structure and abundance at the genus level revealed overall differences between the reference Zlato Pole wetland and the seasonally flooded areas (C and P). It differentiated the non-flooded samples, and the flooded samples, with the lowest level of dissimilarity between non-flooded zones of the paddy fields despite their spatial distance (Figure 2B). The bacterial communities in the flooded rice paddies were positioned in a separate cluster. This is often observed in irrigated rice fields [35], and could explain the considerable number of similar bacterial sequences obtained in the present study.

Figure 2. NMDS ordination of Bray–Curtis dissimilarities based upon differences in microbial community structure (A) and hierarchal cluster analysis at the species level (B).

Note: *nat, reference natural wetland Zlato Pole; art, seasonally flooded rice paddies

The analysis of the microbiomes in the flooded and non-flooded soils confirmed the hypothesis of Foulquier et al. [36] that the taxonomic composition is highly affected by the hydrological regime. At the time of sampling, the seasonally irrigated soils harbored bacterial communities distinct from those in the permanently flooded wetland (ZP), as well as in the dry soils. This is consistent with other findings that hydrology is one of the main drivers that model the structure of the sediment bacterial populations [37]. Further, the bacterial composition of the sediments and adjacent dry soil consortiums in our study demonstrated that the time (5 months) of sediment submersion is insufficient for convergence of communities between permanently and seasonally flooded soils, despite the notable similarity in the environmental conditions. This contradicts the Rees et al. [38] report that microbial communities in sediments subjected to tow month drought, could fully return to their initial structure one month after rewetting. It is possible for the dry-wet cycle rotation in the paddy fields to act as a physiological stressor for the bacterial communities, due to periodical oxygen depletion and alteration in water potential within the sediments [36].

The sediment samples presented lower concentrations of OM, TN and nitrogen forms than those from the adjacent dry areas. The OM content in non-flooded soils stands out as a key factor shaping the bacterial communities, in contrast to water content in the inundated soils. OM content and soil types have been previously described as major driving factors of microbial community structure in terrestrial ecosystems, ensuring the substrate availability for microorganisms [36,39]. In addition, pH varied between the soil types and had a direct effect on the microbial community structure. The pH proved one of the parameters that could explain the observed higher variation in the community among soil types than study sites, as reported in a number of previous studies [40–42]. The seasonally inundated paddy fields were characterized by an accumulation of available phosphorus, as a result of fertilization [43], and harbored higher bacterial diversity than the ZP1 site (Table 2). Hou et al. [44] suggest that paddy-upland reciprocal cultivation is related to rotation of anaerobic and aerobic conditions, affecting the soil C and N cycles. This presumably explains the increase in microbial diversity and suggests that the bacterial diversity in such soils would be higher regardless of the fertilization employed. The results demonstrate that the differentiation in bacterial community compositions of the rice sediments from the adjacent areas and the natural wetland are influenced not only by the flooding but also by the rice production.

Site-specific bacterial community taxonomic composition

The analyzed sequences were distributed among 37 bacterial phyla, with over 97% bacterial sequences classified within only nine of them. The area-specific distribution is presented in the Krona charts (Supplemental Figures S2–S7). The taxonomic composition in the studied sites included some rare and low abundant phyla. The studied communities exhibited taxonomic variation, but were generally dominated by four major phylogenetic groups (Proteobacteria, Actinobacteria, Acidobacteria and Bacteroidetes), accounting for 84.5% of the total number of obtained sequences. The dominant groups are presented in Figure 3A (relative abundance >1%). The bacterial communities in the rice paddies sediments were more evenly distributed, whereas the Zlato Pole sediment was the most biased. Proteobacteria was the most diverse phylum across samples, encompassing 34.2% (P1) to 67.7% (ZP1) of all bacterial sequences, followed by Actinobacteria (7.8%-26.3%). Division Proteobacteria has been frequently reported as a dominant in sediments or soils. Proteobacteria are actively involved in the metabolic processes related to global carbon, nitrogen and sulfur cycling in natural and artificial wetlands [10,11,42,45]. The presence of Actinobacteria was also not surprising, as they are reportedly one of the dominant phyla in a number of soils and wetlands [31,46,47]. Acidobacteria (2%-16.9%) was the third dominant phylum, equally abundant in the paddy field sediments. It is considered a sensitive indicator for environmental changes in response to agricultural activities [48], however there were no specific distribution patterns in the distribution of Acidobacteria subgroups. Further, the Pearson correlation analysis demonstrated no significant correlations with environmental factors (p > 0.05) (Supplemental Table S2). The results were inconsistent with the previous report that members of Acidobacteria are pH sensitive and thrive in acidic environments [42]. In contrast, in the ZP1 sample, Firmicutes (9.4%) and Bacteroidetes (8.3%) dominated over Acidobacteria. The abundance of Chloroflexi and Gemmatimonadetes was also significantly higher (p ≤ 0.05) in the sediments from the rice paddies than in the ZP area. Verrucomicrobia and Nitrospirae were relatively evenly distributed among samples. Other universally represented minor taxa included Cyanobacteria, Chlorobi, Saccharibacteria, Planctomycetes and Tectomicrobia, whereas Deinococcus-Thermus phylum was detected primarily in ZP1 sediments. They are typically found in natural and constructed wetlands [10,31,45], paddy fields [35] and other agricultural soils [49].

Figure 3. Abundances at Phylum level (A) and Class level (B) based on MiSeq Illumina sequencing of the V3–V5 region of the 16S rRNA gene in the wetland samples.

Figure 3B shows the relative abundance in the percentage of bacterial sequences identified at the class level (relative abundance > 1%). In total, 117 bacterial classes were identified, and the abundances of 67 of them were shared between all samples. None of the classes dominated in all communities. Among Proteobacteria phylum, the most abundant classes were Alphaproteobacteria (21.4%), followed by Gammaproteobacteria (12.5%), Betaproteobacteria (6.8%) and Deltaproteobacteria (4%). Alphaproteobacteria and Betaproteobacteria had a more significant contribution to the rice paddies, which is in agreement with previous findings [50]. These taxa are ubiquitous in freshwater sediments, usually correlated with pH and nutrients [13]. Other classes, such as Deltaproteobacteria and Gammaproteobacteria, were represented by a more limited number of sequences compared to other paddy field studies [42]. It has been previously reported that Betaproteobacteria prefer relatively oxygen-rich conditions [32], suggesting that the Tsalapitsa and Plovdiv paddy sediments represented oxic conditions for the bacterial consortium. Moreover, this is supported by the extremely low richness of the Verrucomicrobia phylum than generally observed [13,31], which prefer relatively anoxic conditions [51], and more nitrogen-rich content [52].

By contrast, the microbial consortium in Zlato Pole sediments was characterized by an extremely high abundance of Gammaproteobacteria (45.1%) and Bacilli (Firmicutes). The higher richness of Firmicutes compared to the paddy sediments (p < 0.05) is an indicator of more extreme conditions in the reference wetland. It is consistent with previous observations, and could be explained by their adaptation to thrive in low-nutrient environments through spore-formation [29]. Numerous studies have indicated that Gammaproteobacteria usually thrive in organic-rich environments, such as sites contaminated with agricultural or organic pollution [11,53,54]. They are also involved in reduction processes in anaerobic sediments [10]. Actinobacteria, Acidimicrobia, Acidobacteria subgroup 6, Sphingobacteria and Thermoleophilia, were universally represented. The soils obtained from the non-flooded areas of Zlato Pole wetland, Tsalapitsa and Plovdiv rice paddies displayed higher abundance of Actinobacteria (18.5%, 22%, and 15% respectively) than the samples obtained from the sediments, especially those from ZP1 (6%). Bacilli were ten times more abundant in ZP1 (8.3%) than in the C1 and P1 sediments (0.8%). In addition, some of the minor classes as Deinococci, Flavobacteria and Verrucomicrobiae were presented primarily in the ZP1 community. Others, including Solibacteres, Holophagae, Acidobacteria, Ignavibacteria, etc. were found mainly in the paddy sediments (Figure 4A).

Figure 4. Heatmap analysis of bacterial community at class (total abundance) level (A) and genus (relative abundance >1%) level (B).

Note: *UB, uncultured bacterium.

Of the 885 genera identified in the wetland soils, 357 (40.3%) were detected in all samples. Relative abundance >1% was established for only 51 (5.8%) genera, and the majority of identified bacteria were affiliated to rare taxa. Sphingomonas, Nocardioides, Skermanella, Altererythrobacter, Solirubrobacter, Blastococcus and Nitrospira were detected in high numbers in all samples and were similar for both flooded and non-flooded soils (Figure 4B). The dominant bacterial consortium in the Zlato Pole sediments (ZP1) differed significantly and was strongly dominated by Acinetobacter (53.7%), Sphingomonas (7.4%), Exiguobacterium (6.7%), Cloacibacterium (3%), Massilia (1.8%), Microvirga (1.7%) and Pontibacter (1.4%). Acinetobacter, Exiguobacterium, Cloacibacterium and Pontibacter were detected only in the ZP1 sample. The genus Acinetobacter is able to transform nitrogen by heterotrophic nitrification and aerobic denitrification [55]. It is usually positively correlated with nitrate and total nitrogen content [56]. The high Gammaproteobacteria abundance in ZP sediments, despite the lower nutrient concentration and the absence of correlation with the environmental parameters, is evidence of constant anthropogenic pressure around the area. This suggestion is supported by the presence of Deinococcus-Thermus phylum, recognized as one of the most extremophilic bacterial taxa [57]. Other genera identified only in ZP1 include the sanitary-state indicators Enterobacter, Salmonella and Raoultella. The Tsalapitsa and Plovdiv rice paddies sediment communities (C1 and P1) ware characterized by higher diversity compared to ZP1, including microorganisms actively involved in the nitrogen cycle [10], e.g. genera such as Sphingomonas, Gemmatimonas, Bradyrhizobium, Solibacter, Bryobacter, etc., without the presence of a highly pronounced dominant genus (Figure 4B). The non-flooded soils (C2, P2 and ZP2) had a highly similar bacterial consortium composition. The main differences with the sediment communities were not in the dominant complex, but in the higher abundance of area-specific genera such as Micromonospora, Roseiflexus and Acidibacter. Overall, the dominant complex composition showed similar clustering of the soil communities (Figure 4B) to that based on the total identified taxa (Figure 2B).

Redundancy analysis

Redundancy analysis (RDA) was applied to differentiate the beta-diversity variation based on the set of physicochemical parameters. It resulted in a model that could explain 70.79% of the total variation of the dominant genera in the bacterial community (Figure 5). RDA1 was strongly correlated to C1, P1 and ZP1 sediment communities, and RDA2 was associated with the communities from the non-flooded areas (C2, P2 and ZP2). The bacterial communities in C1 and P1 sediments were associated with soil water content and AP concentration, leading to differences in their taxonomic composition from those in the reference wetland, where, by contrast, the community indicated no relation with the analyzed environmental parameters. In contrast, those in the C2, P2 and ZP2 sites were associated with NO3-N, TN and OM. The communities in the areas adjacent to the non-flooded areas at the C2 and P2 sites, and especially at the ZP1 area, were strongly related to OM and TN concentration. The correlation analysis showed that AP was one of the main driving factors in the bacterial communities at the genus level (Supplemental Table S2). Of the dominant genera, Briobacter, Solibacter, Flavisolibacter and the members of family Anaerolineaceae were positively correlated with AP (p ≤ 0.05).

Figure 5. Relationship between the abundance of dominant genera in the bacterial communities and environmental factors in natural and seasonally flooded wetlands.

The focus in this study was on near-surface soils, since they are more responsive to land management [58]. We used similar experimental conditions for taxon richness estimation. It is critical for the assessment of the effect of the environmental and spatial variation on the soil bacterial species richness within a given region [59]. The bacterial community structures determined by MiSeq Illumina sequencing showed significant differences between the permanently inundated sediments of the Zlato Pole wetland and the seasonally flooded sediments of the protected area “Orizishta Tsalapitsa”. Multivariate analyses demonstrated that the bacterial communities were separated into three groups, suggesting that wetland type, as well as soil type, are the major factors in determining the bacterial community in Maritsa River Basin wetlands. The dominant complex was clearly linked to the NH₄-N, TN and OM content, as reported in a previous study [60]. The soil microbiomes exhibited distinct clustering based on soil type, and land management practice, with paddy field sediment samples clustering together and separately from the adjacent dry soils, and reference wetland sediments, suggesting evidence of a cultivation-specific microbiome, as observed by Lopes et al. [61]. However, since only 70.79% of the total variation of the dominant genera was explained by the RDA, we suggest that there are other factors affecting the bacterial composition in the studied sites, too. Such parameters may include the food web interactions [62], dispersal limitations [63], as well as “top-down” control mechanisms [64].

Conclusions

The present study provides the first detailed analysis of the bacterial diversity in two different wetlands in the Maritsa River Basin, taking into account the effects of the hydrological regime and rice cropping. The bacterial community structures and their differences between natural and seasonally flooded wetlands were clearly distinguished by MiSeq Illumina sequencing. The results revealed a significant difference in the bacterial community structure between the permanently inundated sediments of the Zlato Pole wetland and the seasonally flooded sediments of the protected area “Orizishta Tsalapitsa”. Multivariate analyses suggested that the major factors in determining the bacterial community in Maritsa River Basin wetlands were the type of wetland and the type of soil. The dominant complex was linked to the ammonium nitrogen, total nitrogen and organic matter content. The soil microbiomes clustered based on soil type and land management practice. These observations are important for better understanding the microbial ecology of wetlands, especially considering the water regime as a main driver for the bacterial diversity. The obtained results confirmed the relative importance of environmental parameters on microbial community structure. Nevertheless, further and more detailed study is needed to detect the spatiotemporal patterns of the bacterial communities found in permanently inundated sediments and those found in seasonally dewatered locations.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was funded by the Department for Scientific Research of the Plovdiv University Paisii Hilendarski (grant number FP17-BF-001), and was supported by the National Scientific Program for Young Scientists and Postdoctoral Fellows, Bulgarian Ministry of Education and Science.