Phosphorus forms in the sediment of seagrass meadows affected mainly by fungi rather than bacteria: a preliminary study based on ³¹P-NMR and high-throughput sequencing

Xincun Bay is located in the southeast of Hainan Island, with only one narrow channel connected to the South China Sea. The main driving force behind the water flow is the irregular diurnal tide. The area of the bay is 22.5 km², with an average water depth of 4.2 m and a tidal range of 0.7 m; more than 3.2 km² of the total bay area is allocated for aquaculture. In 2002, seagrasses covered 2.0 km² of sand-mud beaches in Xincun Bay (Huang et al. 2006), but since then the meadows have been declining due to high anthropogenic pressure (Yang & Yang 2009). Sediment samples were collected from monospecific meadows of Enhalus acoroides, Thalassia hemprichii, Cymodocea rotundata, Halophila ovalis and the bare tidal flat in the intertidal zone of Xincun Bay in March 2018 (Fig. 1). Three plots were selected in each treatment. Five replicates of surface sediment samples were randomly collected in each species plot with a small acrylic plastic corer that penetrated to a depth of 5 cm. After the collection, sediment samples were placed on dry ice and immediately transported to the laboratory within 3 h. Some of them were stored at −20°C for 16S and ITS rRNA gene sequencing and the rest was dried by lyophilization for ³¹P-NMR analysis.

Figure 1

A 3.0 g sample of lyophilized sediment was extracted in 30 ml mixture of 0.25 mM NaOH and 0.05 mM Na₂EDTA for 16 h while shaking and then centrifugated (10 min, ~ 10 000 × g). The supernatants were frozen and freeze-dried. The freeze-dried material in the amount of 0.4 g for each sample was re-dissolved in 0.6 ml D₂O and 0.1 ml of 10 mM NaOH, centrifuged (25 min, 12 000 × g) and transferred to a 5 mm NMR tube. The ³¹P-NMR spectra were obtained using the AVANCE IIITM 500 MHz spectrometer (Bruker) at 202.47 MHz. A pulse of 45°, an acquisition time of 0.68 s, a pulse delay of 4.32, a temperature of 20°C, and ~10 000 scans were used. Solution ³¹P-NMR spectra were processed with 5 Hz line-broadening on MestReNova Software (version 12.0.0-20080). To obtain peak areas, different P compounds were identified by their chemical shifts by standardizing the orthophosphate peak to 6.0 ppm. Peaks in the raw spectrum with a signal-to-noise ratio exceeding 7 were fitted using the deconvolution method in MestReNova software. The percentage of individual P compounds was calculated as peak areas relative to the sum of all P peak areas based on previous reports (Shinohara et al. 2012; Baldwin 2013; Bramha et al. 2014).

The steps described below were conducted by Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). Microbial DNA was extracted from sediment samples using the E.Z.N.A.® Soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s protocols. The final DNA concentration and purification were determined by a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA), while DNA quality was checked using 1% agarose gel electrophoresis.

The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with primers 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’) by a thermocycler PCR system (GeneAmp 9700, ABI, USA). PCR reactions for each sample were performed in triplicate in 20 μl mixtures containing 4 μl of 5 × FastPfu Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of each primer (5 μM), 0.4 μl of FastPfu Polymerase and 10 ng of template DNA. The hypervariable regions of the fungal ITS gene were amplified with primers ITS1F (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS2R (5’-GCTGCGTTCTTCATCG ATGC-3’) by the same thermocycler PCR system as above. PCR reactions for each sample were performed in triplicate in 20 μl mixtures containing 2 μl of 10 × Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of each primer (5 μM), 0.2 μl of rTaq Polymerase and 10 ng of template DNA. PCR reactions were conducted using the following program: 3 min of denaturation at 95°C, 27 cycles of 30 sat 95°C, 30 s for annealing at 55°C, 45 s for elongation at 72°C, and a final extension at 72°C for 10 min. The resulting PCR products were extracted from a 2% agarose gel and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, USA) according to the standard manual. Purified amplicons were pooled in equimolar ratios and paired-end sequenced (2 × 300) on the Illumina MiSeq platform (Illumina, San Diego, USA).

Processing of the sequencing data: the steps described below were also conducted/operated by Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). Raw fastq files were demultiplexed, quality filtered by Trimmomatic and merged by FLASH based on the following criteria: (i) reads were truncated at any site that received an average quality score < 20 over a 50 bp sliding window; (ii) primers were exactly matched, allowing for a 2-nucleotide mismatch, and reads containing ambiguous bases were removed; (iii) sequences with overlaps longer than 10 bp were merged according to their overlapping sequence. Operational taxonomic units (OTUs) were clustered with a 97% similarity threshold using UPARSE (version7.1 http://drive5.com/uparse/) and chimeric sequences were identified and removed using UCHIME. The taxonomy of each 16S and ITSrRNA gene sequence was analyzed by the RDP classifier algorithm (http://rdp.cme.msu.edu/) against the 16S rRNA database (Silva 132/16S bacteria) and the ITS rRNA database (UNITE 8.0/its fungi) using a confidence level of 70% (Liu et al. 2019).

Statistical analyses of all parameters were performed using the IBM SPSS statistical software package, version 19 (IBM Corporation, New York, USA). Data from each treatment were analyzed using one-way analysis of variance (ANOVA), and Duncan’s multiple range tests (p < 0.05) were performed for multiple comparisons. Non-metric multidimensional scaling (NMDS) plots were prepared in R with the vegan package (version 3.0.2). Redundancy analysis (RDA) was performed using CANOCO (version 5.0).

The solution ³¹P-NMR spectra revealed six forms of the P compound in the NaOH-EDTA extraction of sediments from various seagrass meadows, where orthophosphate (Ortho-P), pyrophosphate (Pyro-P) and polyphosphate (Poly-P) belong to inorganic P (Pi), orthophosphate monoesters (Mono-P), orthophosphate diesters (Di-P) and phosphonates (Phon-P) belong to organic P (Po) (Fig. 2), based on previously reported chemical shifts (Paytan et al. 2003; Cade-Menun 2015). Pi accounted for 61.36% to 93.40% of the TP extracted, much more than that of Po– from 6.60% to 38.64%. In general, Ortho-P accounted for the highest proportion and was followed by Mono-P, while the remaining forms contributed less in our study.

Figure 2

The variety of P compounds varied in the sediments of different seagrass meadows (Fig. 3). In E. acoroides, all six P forms were present, whereas in T. hemprichii, C. rotundata and H. ovalis – two Pi and two Po compounds were found, in the bare tidal flat – all three Pi forms were found, but only one Po compound. The proportion of P forms in NaOH-EDTA extracts of sediments from different seagrass meadows was also very varied (Fig. 3). The E. acoroides meadows were characterized by the lowest proportion of Pi, while the content of Po was the highest of all samples – up to 38.64%, with Phon-P accounting for 19.05%. The P form and proportion were similar in T. hemprichii and C. rotundata. There was no Di-P in H. ovalis, but a small amount of Phon-P was found. In the sediment of the bare tidal flat, the proportion of Pi was the highest compared to the sediments with seagrasses, up to 93.40%, especially of Poly-P, which accounted for 10.78%. On the other hand, Mono-P accounted for only 6.6% and was the only Po form detected by ³¹P-NMR spectroscopy in the sediment of the bare tidal flat. In general, seagrass sediments contained more forms and a higher proportion of Po compared to the bare tidal flat sediment, especially E. acoroides sediments.

Figure 3

The analysis of 16S and ITS rRNA gene sequences in 15 samples provided a total of 463 575 and 936 900 valid sequences clustered into 2874 OTUs and 1283 OTUs at 97% sequence similarity, respectively. They were assigned to 53 phyla and 514 genera in the case of bacteria, and seven phyla and 356 genera in the case of fungi. In the case of bacteria, there were 1622 core OTUs in five groups; 51, 5, 20, 6 and 4 unique OTUs were identified in the sediments from E. acoroides, T. hemprichii, C. rotundata, H. ovalis, and the bare tidal flat, respectively. In the case of fungi, there were 116 core OTUs in five groups; 197, 120, 151, 107 and 159 unique OTUs were identified in sediments from E. acoroides, T. hemprichii, C. rotundata, H. ovalis and the bare tidal flat, respectively. The ratio of the core OTUs of bacteria was 56.44%, which was much higher than that of fungi (9.04%).

Differences in the diversity of bacterial and fungal communities among all samples are presented in Table 1. The richness indices, Ace and Chao, showed significantly the lowest value for the bacterial community in E. acoroides sediments (2147, 2142), followed by the bare tidal flat and H. ovalis (2392, 2408), while C. rotundata (2424, 2414) and T. hemprichii (2426, 2428) showed the highest richness values. There was no significant difference in the bacterial Shannon and Simpson indices among the samples. With regard to the fungal community diversity in the sediment, the highest values of the Ace and Chao indices were recorded for E. acoroides (303, 303), followed by H. ovalis (256, 259) and the bare tidal flat (246, 246). A significant difference was found in the Shannon index and the Simpson index for fungi. The lowest diversity of the fungal community was recorded for T. hemprichii (2.82, 0.188), followed by H. ovalis (3.87, 0.071) and E. acoroides (3.91, 0.064), while C. rotundata (4.12, 0.040) and the bare tidal flat (4.37, 0.029) showed higher diversity. NMDS plots based on the confidence intervals (p < 0.05) showed no overlap among different sites, which proves that bacterial (stress = 0.105) and fungal communities (stress = 0.087) in the sediment at five sites were different, especially between E. acoroides and other sites of bacterial communities, between the bare tidal flat and other sites of fungal communities (Fig. 4). In other words, E. acoroides can clearly affect sediment bacterial communities, just as seagrasses can affect sediment fungal communities.

Figure 4

The relative abundance of bacterial (> 1% of the total reads) and fungal phyla under different treatments are shown in Figure 5. The dominant bacterial phyla were Proteobacteria (from 29.11% to 38.62%), Bacteroidetes (15.72% to 22.16%), Firmicutes (8.91% to 13.53%), Chloroflexi (8.21% to 11.92%), Actinobacteria (4.65% to 8.96%), Acidobacteria (3.02% to 3.33%), Latescibacteria (1.86% to 2.44%), Cyanobacteria (0.65% to 2.75%), Verrucomicrobia (1.03% to 2.34%) and Planctomycetes (1.08% to 1.92%; Fig. 5A). According to the Duncan test, the relative abundance of Bacteroidetes and Firmicutes was significantly higher in the bare tidal flat than in the seagrasses, and there was no significant difference between different seagrasses. The dominant fungal phyla were Ascomycota (32.56% to 90.51%), Basidiomycota (4.57% to 65.40%) and Rozellomycota (0.14% to 1.07%; Fig. 5B). Interestingly, the relative abundance of Ascomycota was significantly the highest, whereas that of Basidiomycota was significantly the lowest in the bare tidal flat. On the other hand, the relative abundance of Ascomycota was significantly the lowest, while that of Basidiomycota was significantly the highest in the sediment of H. ovalis.

Figure 5

Previous studies explored and screened P-cycling-related bacteria and fungi (Khan et al. 2007; Sharma et al. 2013; Li et al. 2019). At the genus level, only five bacteria and 15 fungi were screened, and the average percentage was respectively 0.29% of the total bacteria and 26.20% of the total fungi. Figure 6 shows the relative abundance of the top five genus-level P-cycling-related bacteria and fungi, which on average accounted for 100% and 97.97% of the total relative abundance of the P-cycling-related bacteria and fungi, respectively. The relative abundance of P-cycling-related bacteria, ranging from 0.12% to 0.60%, was very low and there was no significant difference between the treatments (Fig. 6A). The top five genus-level P-cycling-related bacteria were Bacillus, Streptomyces, Staphylococcus, Pseudomonas, and Paenibacillus, which belong mainly to the Bacilli class and only Streptomyces belongs to Actinobacteria. Compared with bacteria, the relative abundance of P-cycling-related fungi (15.98% to 39.31%) was higher, with H. ovalis and the bare tidal flat having the highest abundance, followed by E. acoroides and C. rotundata, T. hemprichii with the lowest abundance (Fig. 6B). The top five genus-level P-cycling-related fungi were Cladosporium, Penicillium, Aspergillus, Trichoderma and Candida, belonging to different classes but the same phylum Ascomycota. The dominant fungi were Cladosporium and, according to the Duncan test, the relative abundance of H. ovalis was significantly higher compared to other treatments. The lowest relative abundance of Penicillium was observed in T. hemprichii, while the highest in H. ovalis, and the difference between them was statistically significant. The highest relative abundance of Aspergillus was observed in the bare tidal flat and was significantly higher compared to other treatments. The highest relative abundance of Trichoderma was recorded in E. acoroides, while the lowest in T. hemprichii, and the difference between them was significant. In addition, the highest relative abundance of Candida was recorded in E. acoroides, and the lowest in the bare tidal flat, and the difference between them was significant.

Figure 6

Redundancy analysis and a Monte Carlo permutation test were performed to determine the relationships between the top five genus-level P-cycling-related microbes and the sediment P forms. Visualized through RDA ordination (Fig. 7), all of the P form variables explained statistically 38.31% of the bacterial variance, with axis 1 explaining 38.23% of the variance and axis 2 explaining 0.08% (Fig. 7A). One of the issues emerging from these findings is that different forms of P have different effects on the abundance of P-cycling-related bacteria. Of the P form variables considered in this study, the Pyro-P negatively correlated with Bacillus and significantly explained the variations of P-cycling-related bacteria (p < 0.05), indicating that this was the principal factor affecting the P-cycling-related bacterial communities. All of the P fraction variables explained 73.63% of the fungal variance, with axis 1 explaining 61.90% of the variance and axis 2 explaining 11.73% (Fig. 7B). Mono-P and Di-P were suggested to significantly explain the variability in P-cycling-related fungi composition (p < 0.05). Mono-P, the main factor determining the distribution of P-cycling-related fungi, was negatively correlated with Aspergillus and Trichoderma. Di-P, another important factor, was positively correlated with Candida and negatively correlated with Penicillium and Cladosporium. By comparison, the other factors had no obvious effect on P-cycling-related fungal communities. The main finding of this research is presented in Figure 8.

Figure 7

Figure 8

To our knowledge, little data has been published so far on P forms using ³¹P-NMR in seagrass meadows. The results from ³¹P-NMR spectroscopy showed that Ortho-P was the most dominant P compound in all our samples. The proportion of Ortho-P, accounting for 58.61% to 81.70% in Xincun Bay sediments, was similar to other marine and estuary sediments (Watson et al. 2018; Prüter et al. 2019). The Poly-P signal was not observed by the solution-state NMR spectra in the seagrass meadows, except for E. acoroides. Poly-P was considered a labile form of P with a half-life of 0.8 to 2.0 years and could be transformed to Pyro-P or tri-metal-phosphate catalyzed by various metal ions, which consequently was rarely detectable by the solution-state ³¹P-NMR in sediments (Liu et al. 2009; Ding et al. 2010; Xie et al. 2019). Pyro-P was observed in all samples and its highest percentage was found in the surface sediments of E. acoroides, which could indicate high activity of soil fungi in the transformation of Mono-P and Di-P into Pyro-P (Koukol et al. 2008).

Mono-P was the dominant Po form in these sediments; it is derived from plants and was detected in all samples. Degradation and mineralization of Mono-P is an important source of net primary productivity (Yuan et al. 2015). The usual Mono-P species to be determined are myo-inositol hexakisphosphate and its stereoisomers, such as D-chiro-inositol hexakisphosphate neoinositol hexakisphosphate and scyllo-inositol hexakisphosphate (Turner et al. 2003). Some of the Mono-P species may be unidentified Di-P degradation products (Schneider et al. 2016). A small portion of Di-P was detected in E. acoroides, T. hemprichii and C. rotundata. Being a type of unstable Po compound, Di-P is usually degraded within a few days and mineralized faster than Mono-P, while playing a key role in the P cycle (Dell'Anno & Danovaro 2005; Reitzel et al. 2006). In addition, Di-P can be degraded to Mono-P during sample preparation and extraction (Turner et al. 2003), so caution should be taken in the interpretation as these concentrations do not necessarily reflect the actual content in the sediments. The observed absence or low proportion of diesters, ranging from 0% to 4.98%, in the sediments could be attributed to easy degradation in nature and the extraction time of 16 h in this study. Phon-P in the present study appeared to be greatly variable, ranging from 0% to 19.05% in different samples, which is consistent with many other reports (Benitez-Nelson et al. 2004; Sannigrahi & Ingall 2005). The highest Phon-P proportion of 19.05% in E. acoroides sediments is probably explained by the interception of organophosphorus herbicides from farmland and insecticides from the nearby aquaculture. It is not possible for the secretion and degradation of seagrasses to produce such a high Phon-P concentration naturally.

The present study has shown species-specific differences in P forms in seagrass meadows, especially in E. acoroides. Yuan et al. (2015) reported that the proportion of organic P in total P in the sediments varied in different plants. On the one hand, according to the existing studies (Duarte & Chiscano 1999) and our research, the biomass and production of E. acoroides are much higher compared to other seagrasses, so litter production of E. acoroides is also higher. On the other hand, E. acoroides, the largest (over 1 m high) seagrass species in Xincun Bay, can capture more different floating organic matter from water into sediments by slowing down the water flow (Komatsu et al. 2004). The vertical carbon flux can be estimated at 557.6 g C m⁻² yr⁻¹ in E. acoroides meadows (Wahyudi et al. 2016). The higher concentration and greater variety of organic matter in E. acoroides meadows result in higher organic P content and more abundant P forms compared to other seagrass meadows.

Tropical seagrasses, generally P-limited owing to the strong P fixation capacity of carbonate-rich sediments, can surprisingly form densely vegetated meadows in such low-nutrient environments (Brodersen et al. 2017). The high P mobilization in tropical seagrass sediments could potentially be further supported by microorganisms (Vazquez et al. 2000). Microorganisms have been increasingly recognized as pivotal players in seagrass ecology (Ugarelli et al. 2017; Brodersen et al. 2018). Previous studies focused mainly on the effect of microorganisms on carbon, nitrogen and sulfur cycles, but research on the effects on the P cycle is limited. Jiang et al. (2015) reported that T. hemprichii hosted the highest abundance of bacteria, followed by E. acoroides and C. rotundata, but in this study H. ovalis was characterized by the highest abundance of bacteria. Similarly to the findings of Hurtado-McCormick et al. (2019) regarding the Zostera muelleri seagrass microbiome in the sediment where Flavobacteria, Chronaxies, and Desulfobacteraceae dominated the bacteria communities at the order level, we found that Flavobacteriaceae, Clostridiales and Desulfobacteraceae were the main orders and they were present in similar numbers. The finding presented in this study that Ascomycota and Basidiomycota dominated in the E. acoroides sediment thus far supports the idea of Wainwright et al. (2019) at the fungal phyla level. Contrary to the previous studies, it is somewhat surprising that high abundance of Tremellomycetes and Cystobasidiomycetes was detected at the class level, 45.97% and 16.38% in the T. hemprichii sediment respectively, and very low abundance, ranging from 0% to 1.24%, in the bare tidal flat sediment. Physicochemical parameters of the sediments are significantly modified by seagrasses. The bacterial abundance and biomass in the rhizosphere is about twice as high as in the nearby non-vegetated sediments (Jiang et al. 2015; Ugarelli et al. 2017). The seagrass rhizosphere microbiome shows statistically significant differences in the composition compared to adjacent vegetated sediment communities (Cucio et al. 2016; Ettinger et al. 2017). This case study confirms the meaningful difference in the microbial community between the bare tidal flat and seagrass meadows, especially in the case of the fungal community between treatments.

In the natural environment, numerous bacteria and fungi in the soil and sediment, playing an important role in the P cycle, are effective in releasing P from total P through acidification, secretion of organic acids or protons, chelation exchange reaction, mineralization of Po by acid phosphatases, phytases, phosphonates and C-P lyases (Sharma et al. 2013; Dipta et al. 2019; Y. Li et al. 2019). In seagrass sediments, the dominant genus of P-cycling-related bacteria was Bacillus, followed by Streptomyces and Staphylococcus (Fig. 6). These results are similar to those obtained in the study by Li et al. (2019) on the sediments from Lake Sancha, who also found that the most phosphate-solubilizing bacteria strains belonged to the genus Bacillus, followed by Paenibacillus. Bacteria with a high potential for solubilizing phosphate from the rhizosphere of marine plants Avicennia belong mainly to Bacillus, Pseudomonas and Acinetobacter (Teymouri et al. 2016). Of all microorganisms, it was found that fungi have a greater capacity to solubilize insoluble phosphate than bacteria (Sharma et al. 2013; Dipta et al. 2019). In this study, the proportion of P-cycling-related fungi is much higher than that of P-cycling-related bacteria, ranging from 26.20% to 0.29%. Some soil fungi species are well known for their strong potential for P-solubilization or mobilization of Po in soils and sediments, such as Penicillium and Aspergillus (Mercl et al. 2020). Penicillium species are considered a key group of fungi in the P cycling due to their ability to solubilize inorganic P through the release of organic anions (Richardson et al. 2011). Consistent with the aforementioned literature, this study also found that Cladosporium, Penicillium and Aspergillus were the main fungi at the genus level in the seagrass meadows.

Po is an important component of the P pool, but it cannot be used directly by plants and microorganisms. The degradability of Po compounds depends mainly on physicochemical and biochemical properties of their molecules, nucleic acids and phospholipids. Sugar phosphates decompose easily, whereas phytic acid, polyphosphates and phosphonates decompose relatively slowly (Dipta et al. 2019). Previous studies on P-cycling-related microorganisms focused mainly on Pi-cycling-related bacteria and only a few studies investigated Po-cycling-related bacteria and fungi that can excrete phosphatase to degrade organic phosphorus, such as phytates, phosphomonoesters and phosphotriesterases, thereby enhancing the available soil phosphorus level (Khan et al. 2007; Sun et al. 2017). Different microorganisms secrete different enzymes. According to the summary by Dipta et al. (2019), Bacillus, Aspergillus and Trichoderma discussed in this study can secrete alkaline phosphatase, acid phosphatase and phytase; Pseudomonas and Penicillium can secrete alkaline phosphatase and acid phosphatase; Streptomyces can secrete phytase, while Aspergillus and Trichoderma can secrete phytase. Phytate (myo-inositol hexakisphosphate) was quantitatively the most important Mono-P. Therefore, Mono-P was negatively correlated with Aspergillus and Trichoderma in this study.