Biochar amendments in a California salt marsh restoration reduced denitrification and supported distinct microbial community functions

Accelerated sea level rise, combined with human-induced changes to hydrology and sediment transport pathways poses an existential threat to the survival of coastal wetlands over the next century. To support the preservation of wetlands, restoration and enhancement strategies—such as adding sediment to marsh surfaces to boost their resilience to sea level rise—are being actively tested and monitored. This study aimed to assess how biochar derived from Eucalyptus spp. feedstock influences sediment properties, microbial communities, functional potential, and nitrate reduction pathways in a tidal wetland restoration project. Biochar-amended sediments showed reduced denitrification rates, no detectable DNRA activity, and only subtle changes in sediment properties and microbial nitrogen cycling, with the exception of increased carbon content and higher relative abundance of Cyanobacteria and Truepera, a member of the Deinococcota.  While there was some evidence of increased microbial abundance, biochar did not enhance nitrogen reduction as hypothesized, and in some cases appeared to reduce microbial diversity. While biochar additions may benefit the restored wetland by increasing soil organic content to enhance microbial activity, and contributing to carbon sequestration through direct carbon burial, our results highlight the importance of different short-term vs. longer-term impacts to function, and the value of a priori testing of biochar to determine if the type of biochar chosen will chemically react with sediments as intended, before large-scale incorporation into the restoration project.

Methods
Research was conducted at Elkhorn Slough (36.8129°N, -121.7556°), a 12,000-ha estuary tidal estuary in Central California, which was the site of an ambitious 24-ha sediment placement and salt marsh restoration project completed in 2018 (Fig. 1) (termed "Hester Marsh"; Fig. 1). During 2017-8, 180,000 m3 of sediment was used to rebuild a low elevation marsh plain, which was diked and drained in the mid-20th century. Sediments included material from Pajaro River floodplains along with upland soils from a nearby regraded hillside. Key restoration goals were to recover lost ecosystem functions such as carbon storage, nutrient sequestration, and denitrification.

Biochar sediment amendments (10% v/v) were incorporated in 2019 in ten sets of paired plots across the mid-marsh platform. Biochar plots (1 x 1 m) received 5000 cm3 of Eucalyptus spp. biochar mixed into the 5-10 cm sediment layer, capped with 5cm of untreated sediment to prevent washout (Fig. 1). Paired (1 x 1 m) plots were disturbed in a similar manner. Biochar was produced on site by “top-down” pyrolysis of Eucalyptus spp. wood stock, in which the wood material was ignited at the top of the pyre to produce high temperatures that result in an anoxic center as the fire burns downward. Biochar was crushed and sorted, and the <1.5 cm particle diameter was utilized. A nearby natural marsh ("Yampah Marsh") was chosen as a reference marsh due to its spatial proximity to the restoration site.

In October 2021, 160 cm3 sediments samples (four composited 40 cm3 samples from 5-10 cm depth) were collected from biochar-amended (n=10) paired, unamended restoration plots (n=10), and from the nearby natural marsh (n=5; Fig. 1). We also composited fresh biochar and unamended restoration sediments (n=5). Immediately upon return to the lab, natural marsh, and restoration sediments were subsampled for DNA extraction and frozen at -80 °C. Samples were immediately homogenized, and processed for denitrification and DNRA potential activity. The remaining sediment was kept at 4°C prior to physicochemical characterization.

To evaluate potential denitrification and DNRA rates, slurry incubations with 15N tracer were conducted on 15 restoration soils (five unamended restoration samples, five freshly amended biochar samples, five aged biochar amended samples), and five samples collected from the nearby natural marsh. Subsamples of field-wet sediment (2.5 ± 0.5 g) were mixed with 10 mL artificial seawater  purged with argon to remove oxygen, along with glass bead to assist in homogenization. Slurries were mixed in four replicate 12 mL borosilicate Exetainer vials (EA Consumables, LLC, Marlton, NJ), to conduct paired time point measures immediately after injection of the tracer and 18 hours (times T0 and T18, respectively) of incubation for both denitrification and DNRA potential (two for denitrification; two for DNRA). Vials were sealed without headspace and vortexed for 3 minutes. Slurries were pre-incubated in the dark at ambient temperature overnight to remove background 14NOX- and to ensure samples were anoxic. After pre-incubation, vials were injected with 100 µM Ar-purged Na15NO3 solution (98%, Cambridge Isotope Laboratories, Inc., Andover, MA). To halt microbial activity, slurries were injected with 200 µL Ar-purged ZnCl2 solution on ice: half of the tubes were sacrificed at time 0, and half after 18 hours of incubation in the dark. Vials were subsequently stored at 4 °C submerged in artificial seawater until analysis. Dissolved gases (30N2 and 29N2) were measured on a quadrupole MIMS with an inline furnace at the Tropical Research and Education Center, University of Florida.

Potential denitrification rates for each sediment type and treatment were calculated following the isotope pairing technique (Nielsen, 1992). Concentrations of dissolved gases 29N­2 and 30N2 were corrected for sediment and water volumes. Then, differences in concentrations were used to determine production rates of 29N2 and 30N2, or p29 and p30, respectively, normalized by the weights of sediment subsamples in each vial. Rates of denitrification of 15NO3- (D15) and 14NO3- (D14), were determined using  the following equations to calculate D14 and D15:

D15 = p29 + 2p30

D14 = D15[p29/2p30]   

Finally, D14 and D15 were summed for total denitrification rate (Dtotal) for each sediment type and treatment. DNRA potential was measured using the OX/MIMS method. Half of the replicates for DNRA measurements were oxidized with hypobromite iodine solution prior to analysis of 29+30N2 on a MIMS with an inline furnace, to convert NH4+ to N2. . The difference between unoxidized and oxidized samples represented the concentration of 15NH4+. Rates were then calculated as production of 15NH4+ using the same methods described above.

Sediment organic content, percent water content, and the field-moist/oven-dry ratio were determined by drying a 5 cm3 sediment volume for 24 hours at 105 °C, followed by combustion at 550 °C for 5 hours. Sediment pH was measured using a Thermo Scientific Orion A121 meter on vortexed 1:1 (w/w) deionized water to dry sediment slurries. Slurries were produced without pre-drying samples based on the field-moist to oven-dried ratio measured above. Sediment salinity was measured using a YSI Pro Series multimeter with conductivity probe on vortexed and centrifuged 10:1 DI to dry sediment slurries. Measured salinities were adjusted to a 1:1 water:sediment ratio. Soils were extracted with a 2M KCl solution, and extractable NH4-N and NOx--N were measured on a Lachat QC8500 Flow Injection Colorimeter at Stony Brook University's Center for Clean Water Technology. Samples were additionally measured for redox potential (herein, referred to as eH) by placing the probe in sediment to a depth of 5-cm using a benchtop Oakton ORP electrode after allowing readings to stabilize for at least 10 minutes; readings were corrected using standard solution.

Genomic DNA was isolated from 23 samples (10 from restoration sediments, 10 from biochar amended restoration sediments, and 3 from reference marsh) using a Takara Nucleospin Soil DNA extraction kit (Takara Bio USA, Inc., San Jose, CA) following manufacturer’s instructions. Extracted DNA was quantified with a Qubit 3.0 fluorometer. Library preparation and amplicon sequencing were conducted at the University of Minnesota Genomics Center. The amplicon library was generated for V3-V5 region of the 16S rDNA gene following protocols outlined in Gohl et al. (2016).  Amplicon libraries were sequenced using 2 x 300 PE sequencing chemistry on the MiSeq® platform (Illumina). The raw sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) database under accession number PRJNA1263012.

Bioinformatics and statistical analysis were conducted using R v4.3.1. Raw paired-end reads were primer-trimmed using Cutadapt  by calling Python v. 3.1.1 through R. Upon reads’ quality inspection, only the forward reads were used to generate amplicon sequence variants (ASVs) using “dada2” v1.28.0. Chimeras and singletons were filtered and produced amplicon sequence variants (ASVs) taxonomically classified using Silva reference database release 138. 

Data matrices were managed using R package “phyloseq” v1.28.0 and plots were generated using R package “ggplot2” v3.4.2. Amplicon abundance data were rarefied to 500 reads per sample prior to data analysis. Principal Component Analysis (PCA) was carried with edaphic parameters to visualize differences among samples from natural and restored marsh and among control, fresh biochar and aged biochar plots. Non-metric multidimensional scaling (NMDS) analyses coupled with the “envfit” procedure that computes correlations of environmental factors with NMDS axes explored patterns of microbial assemblage composition. PCA and NMDS ordinations were conducted using R package “vegan” v2.6-4. Alpha diversity metrics were calculated to compare diversity in biochar-amended vs. restoration and reference plots as indicators of richness (Chao1), dominance (Berger-Parker), phylogenetic diversity (Faith's PD), and entropy (Shannon). Diversity indices were calculated using the packages "BiodiversityR", "vegan" and "phyloseq". To calculate Faith's PD, an alignment of bacterial sequences was generated and a phylogenetic tree constructed from these using package "phangorn".

To functionally annotate Amplicon Sequence Variants (ASVs), the ASV table was formatted into appropriate input files for Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2). The ASV table, generated from DADA2, was converted into BIOM (.biom) and FASTA (.fasta) file formats using R v4.3.1. The BIOM format retains the sample-by-ASV abundance matrix, while the FASTA file provides nucleotide sequences of ASVs, both of which serve as essential inputs for PICRUSt2. These files were transferred to the Stony Brook Seawulf cluster, where PICRUSt2 was executed in a Bash shell terminal.

PICRUSt2 predicts the functional potential of microbial communities by placing ASVs onto a reference phylogeny and inferring gene families and metabolic functions based on ancestral trait reconstruction. The output includes Enzyme Commission (EC) numbers and KEGG Orthology (KO) annotations, providing insight into the functional capabilities of the sampled microbial communities. The functional prediction outputs from PICRUSt2 (EC and KO abundances per sample) were analyzed in R to explore differences in functional composition across samples. To visualize pathway abundance across samples, the PICRUSt2 output was filtered to retain only pathways associated with decomposition, followed by total sum scaling (TSS) normalization. One sample (biochar amended sediment, site 7; BcharP7_S165_L001_R1_001.fastq.gz, or) was removed due to low quality. The resulting dataset was organized into a long-format table containing sample IDs, pathway names, experimental group assignments (natural marsh, restoration, and biochar amended), and TSS-normalized relative abundances.

Functional beta-diversity was further examined using non-metric multidimensional scaling (NMDS), implemented with the metaMDS() function from the vegan package. Ordination plots were produced in two dimensions, with ellipses illustrating group-level dispersion. Permutational Multivariate Analysis of Variance (PERMANOVA) was used to evaluate whether differences in functional composition were significant among experimental groups. This was performed using the adonis2() function in the vegan package based on Bray-Curtis dissimilarities. Statistical significance was assessed using 999 permutations. Post hoc pairwise comparisons were conducted following a significant global PERMANOVA result using the pairwise.adonis2() function from the pairwiseAdonis package. P-values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) method.

One-way Analysis of Variance (ANOVA) tests was used to compare denitrification potential rates between sample types (unamended restoration, freshly amended restoration, aged biochar amended restoration, and natural marsh sediments), after log10 transforming prior to analysis to meet assumptions for parametric tests. Non-parametric Kruskal-Wallis rank sum tests were conducted on KO/KEGG functions associated with nitrogen metabolism, alpha-diversity metrics, and edaphic parameters measured (pH, salinity, eH, moisture, LOI, NH4+-N, and NOx--N). For edaphic characteristics, a Bonferroni correction was used, resulting in a corrected-α of 0.0063. Significant interactions were followed by using a post hoc Tukey’s HSD Test (R package “multcomp” v1.4.23) for parametric analyses and a Dunn Test (R package “FSA” v0.9.4) for non-parametric analyses.

Correlation matrices were used to visualize relationships between edaphic parameters and denitrification potential (“corrplot” v0.92). Significance of correlation coefficients was determined using the function ‘cor.mtest’ and a confidence interval of 0.95. Resulting p-values (Supplemental Material) were considered significant if they were less than the Bonferroni-corrected α (p < 0.0014) due to the multiple comparisons and possible Type I error.