Shotgun Metagenomics Reveals the Benthic Microbial Community Response to Plastic and Bioplastic in a Coastal Marine Environment

1Microorganisms rapidly colonize and form biofilms on biotic and abiotic surfaces in marine environments. As such, marine microbial communities are commonly distinguished as free-living or particle-associated. A particle-associated lifestyle is thought to provide microorganisms with increased nutritional resources and environmental stability. Detrital aggregates, in particular, are hot spots of microbial diversity that are fundamentally distinct in comparison to free-living communities. Particles are also habitats of enhanced enzyme activity that contribute to the rapid turnover of particulate organic matter and it is clear that particle-associated microorganisms play an important role in biogeochemical processes.

2Plastic pollution has introduced a new surface for microbial colonization and biofilm formation. A recent global assessment of all mass-produced plastics estimated that 4.9 billion tons of plastic waste was discarded in landfills or natural environments. Oceans serve as a sink for said waste with an estimated 4.8–12.7 million metric tons entering the oceans annually. It is clear that marine microorganisms colonize and form biofilms on floating plastic debris and it is also clear that plastic-associated microbial communities are distinct compared to free-living communities.

3Several studies have characterized pelagic microbe-plastic interactions but relatively few studies have characterized microbe-plastic interactions in benthic systems. In the North Pacific Gyre, a previous study reported that microplastics were colonized by diatoms and a separate study reported an increased concentration of Chl a and an increased abundance of nitrogenase genes in microplastic biofilms. By contrast, understanding of benthic microbe-plastic interactions is limited. The transport of plastic debris to benthic systems is facilitated by sinking; debris with a density higher than seawater (1.02 g cm−3) sinks immediately while floating debris loses buoyancy with biofouling. Sinking is therefore an important process and benthic systems are sinks for plastic debris accumulation.

4In coastal regions, benthic microorganisms support primary producers in the overlying water column by the remineralization of organic-rich debris. The remineralization of carbon in anoxic coastal sediments is predominantly carried out via sulfate reduction. The microbial pathway for dissimilatory sulfate reduction involves the reduction of sulfate to sulfite by sulfate adenyltransferases (Sat) and adenylyl-sulfate reductases (AprBA), followed by the reduction of sulfite by sulfite reductases (DsrAB). Climate warming and eutrophication are expanding the distribution of anoxic sediments worldwide and recent studies have clearly shown that sulfate and sulfite reducers are more abundant, diverse, and widespread than previously believed. Yet, the effect plastic pollution may have on sulfate-reducing microorganisms (SRM) is an open question.

5In 2018, the United Nations Environmental Program reported that more than 60 countries approved bans or levies on single-use plastics. In the same year, the European Parliament approved a ban on the top 10 single-use plastics and called for a reduction in single-use plastics. The proposal states “this systemic change and material substitution will also promote bio-based alternatives and an innovative bioeconomy.” In the wake of legislative actions and increased public awareness, the production of bioplastics such as polyhydroxyalkanoate (PHA) is predicted to increase tenfold. If bioplastic is similarly discarded, the subsequent increase in bioplastic waste entering the oceans will introduce yet another surface for microbial colonization; and, while some research has investigated the response of marine organisms to bioplastic, little is known about how microorganisms respond to bioplastic versus petroleum-based plastic.

6The aim of this study was to characterize the microbial communities and individuals that form biofilms on plastic (polyethylene terephthalate; PET) and bioplastic (PHA) in a coastal benthic habitat. For comparative purposes, ceramic pellets were included as a biofilm control. This study was primarily designed to address two questions. Does the introduction of either plastic or bioplastic select for a distinct microbial community compared to a biofilm control, and does this introduction promote a significant shift in taxa or enzyme pools with implications for polymer degradation or biogeochemical cycling? We hypothesized that the PET biofilms would be indistinguishable from the ceramic biofilms whereas the introduction of PHA would select for the growth of a distinct microbial assemblage involved in polymer degradation and biogeochemical cycling.

7The Laguna Madre (Texas, USA) is a large coastal lagoon in the northern Gulf of Mexico (NGOM). The lagoon is divided into two sections: Upper Laguna Madre (ULM) and Lower Laguna Madre (LLM). The microcosms described in this study were deployed adjacent to a dredge material island at 27°32′39.0“N and 97°17’07.7”W in the ULM (Supplementary Figure S1).The ULM is 76 km in length, 6 km in width, has an average depth of 0.8 m, and is separated from the NGOM by the Padre Island National Seashore. The dredge material island itself is remote to the nearest population center (Corpus Christi, TX) and is therefore relatively free of plastic pollution.

8A microcosm was designed to circumvent the challenges of isolating preexisting plastic debris from heterogeneous sediment samples. A microcosm also permitted the inclusion of a ceramic biofilm control. Briefly, 6.0 g of ceramic pellets (Lyman Products, Middletown, CT, USA), 3.0 g of PET pellets (M&G Chemicals, Ettelbruck, Luxembourg), and 3.0 g of PHA pellets (Doctors Foster and Smith, Rhinelander, WI, USA) were deployed in triplicate inside custom-made microbial capsules (MicroCaps; Supplementary Figure S2) at the sediment–water interface for 28 days (May 24, 2016 to June 21, 2016). All pellets were approximately 3–4 mm in diameter and, therefore, the PET and PHA pellets were considered microplastics. MicroCaps utilized 315-μm Nitex mesh to contain the pellets and permit the exchange of water, nutrients, bacteria, and some grazers but limit the entry of larger organisms. The microcosms and triplicate 1-L seawater samples (also from the sediment–water interface) were collected at the end of the exposure. All samples were stored on ice and processed immediately upon return to the laboratory, with the extraction procedure starting within 2 h of sample collection. Substrates were washed three times with 25 ml of 0.22-μm filter-sterilized, site-specific seawater to remove any organisms not part of the biofilm. Seawater (100 ml) was pre-filtered through 315-μm Nitex mesh and the microbial community was collected by membrane filtration on a 0.22-μm polycarbonate filter (Millipore Sigma, Burlington, MA, USA) in triplicate. Environmental parameters from the start and end of the exposure period were collected using a 6,920 V2–2 Multi-Parameter Water Quality Sonde (YSI, Yellow Springs, OH, USA) and varied minimally as follows: water temperature (28.56–29.28°C), salinity (37.78–40.28), pH (8.35–8.47), and DO (8.30–5.68 mg L−1).

9Genomic DNA was isolated in triplicate for each sample type (seawater, ceramic, PET, and PHA) using a modified version of a high-salt and sodium dodecyl sulfate-based method. The only modification from the original procedure was that instead of 5.0 g of soil, DNA was isolated from a quartered 0.22-μm polycarbonate filter (seawater), 6.0 g of ceramic pellets, 3.0 g of PET, or 3.0 g PHA pellets. Due to their higher density, 6.0 g of ceramic was approximately equal in volume to 3.0 g of either plastic type. Following isolation, DNA was quantified (ng μl−1) and assayed for quality (A260/A280 and A260/A230) using a BioPhotometer D30 (Eppendorf, Hamburg, Germany) and the “ds_DNA” methods group with default settings. Final DNA concentrations were verified using a Qubit Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The DNA was stored in the dark at −20°C prior to sequencing.

10Metagenomic library preparation and sequencing were carried out by Molecular Research LP (Shallowater, TX, USA). A total of 12 metagenomes were sequenced: three seawater, three ceramic, three PET, and three PHA. Libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina Inc., San Diego, CA, USA) and diluted to 10.0 pM. The average library size was determined using an Agilent 2,100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed on the Illumina HiSeq 2,500 platform using 2×150 bp paired-end read chemistry. The DNA concentration (ng μl−1) and average size (bp) of the sequencing libraries and the number of sequence reads produced are reported in Supplementary Table S1. Raw sequence reads were submitted to the European Nucleotide Archive and the European Bioinformatics Institute’s (EBI) Metagenomics Pipeline version 3.0, which includes an automated workflow for read processing. Briefly, overlapping reads were merged with SeqPrep, low-quality reads and adapter sequences were trimmed using trimmomatic, and reads less than 100 bp were removed using BioPython. Identification of the 16S rRNA reads was performed using rRNASelector and FragGeneScan was used to find reads containing predicted coding sequences (pCDS) greater than 60 nucleotides in length. The number of sequenced reads ranged from 26,860,030 to 36,522,626, with an average of 30,888,258 reads per sample. Following quality control and merging, the number of remaining reads ranged from 11,273,807 to 22,918,149, with an average of 16,270,447 processed reads per sample.

11Operational taxonomic units (OTUs) were assigned using the QIIME version 1.9 closed-reference OTU picking protocol and the SILVA 128 SSU 97% database with reverse strand matching enabled. Beta-diversity was analyzed using weighted UniFrac values calculated in QIIME. Permutational Multivariate Analysis of Variance (PERMANOVA) was used to test for significant differences between communities using Primer 7 with the PERMANOVA+ package (PRIMER-E Ltd., Plymouth, UK). PERMANOVA was performed using 999 permutations based on the weighted UniFrac from the beta-diversity analysis in QIIME. Due to the low number of unique permutations possible for pairwise tests, Monte Carlo simulations were used to generate p’s for all pairwise comparisons.

12Predicted coding sequences (pCDS) from each of the 12 metagenomes were aligned against a modified database1 of protein sequences from the ESTHER database of alpha/beta hydrolases. Proteins included in the modified database were restricted to depolymerases, esterases, lipases, and cutinases as these were previously implicated in plastic-microbe interactions. To reduce redundancy within the database, CD-HIT version 4.7 was used to cluster protein sequences with default settings. The modified database used here included 10 protein families containing a total of 1,079 protein sequences. Alignments were performed using the command-line version of blastp, as implemented in BLAST+ version 2.6.0, with an expect value (e) cutoff of 10−5, a minimum alignment length of 30 amino acids, and a minimum percent identity of 50%. To visualize the effect of sample type on hydrolase profiles, proportions of positive alignments were normalized to the total number of positive alignments per sample, and a hierarchical cluster analysis based on a Bray-Curtis resemblance matrix was performed. To test for significant differences among profiles, a similarity profile analysis (SIMPROF) using 999 permutations and a significance level of 0.05 was performed in Primer 7. Additionally, the abundance of protein sequences was normalized to the total pCDS and compared between all samples (n = 12).

13Processed reads for each of the PHA samples (n = 3) were co-assembled using MEGAHIT version 1.1.1. Co-assembly was performed using the “meta-large” preset for large and complex communities, and MEGAHIT was called as follows:

14megahit –presets meta-large -r input --min-contig-len 1,000 -o output -t 40.

15To recover the coding sequences of a significantly enriched polyhydroxybutyrate (PHB) depolymerase (identified using the hydrolase database above), the representative PHB depolymerase gene sequence was aligned to the 118,520 co-assembled contigs using tblastn with an e cutoff of 10−5, a minimum alignment length of 100 bp, and a minimum percent identity of 50%. An alignment to NCBI’s RefSeq genomic database was carried out with the same parameters in an effort to include gene sequences from previously described organisms. The relatedness of the 46 PHB depolymerase gene sequences from this study and gene sequences from the top 10 most similar genera from RefSeq was inferred by constructing a maximum-likelihood (ML) tree with IQ-TREE version 1.6.1. Sequences were aligned using the M-Coffee mode of T-Coffee that combines results from eight different aligners and the tree was generated with 1,000 ultrafast bootstraps using the best-fit model (TPM3u + F + R3) as determined by ModelFinder. The phylogenetic tree was annotated with FigTree version 1.4.32. Sequence similarity between the 46 PHB depolymerase gene sequences was determined using trimAl version 1.4.15.

16To characterize individual SRM within the PHA biofilm community, metagenome-assembled genomes (MAGs) were recovered using a previously described method. Briefly, the processed reads were mapped to the co-assembled PHA metagenomic contigs (produced using MEGAHIT above) with BWA version 0.7.15-r1142 using the BWA-MEM algorithm with default parameters. Genomes were recovered with MetaBAT version 2.12.1 using default MetaBat2 settings and a minimum contig size of 2000 bp. The resulting bins were merged, filtered, and refined using CheckM version 1.0.11 and RefineM version 0.0.23, as described previously. Of the 46 bins produced, six MAGs were retained for further exploration based on criteria adopted by Parks et al.: an estimated quality of ≥50 (completeness – 5× contamination), scaffolds resulting in an N50 ≥ 10 kb, containing <100 kb ambiguous bases, and consisting of <1,000 contigs and < 500 scaffolds. Of those six MAGs, the three belonging to SRM as inferred with the “tree_qa” function in CheckM were explored further.

17The three SRM MAGs were analyzed with the Pathosystems Resource Integration Center’s (PATRIC) comprehensive genome analysis service. The PATRIC database contains over 190,000 bacterial genomes and has been increasingly used in environmental studies.The automated PATRIC service includes annotation with RASTtk, prediction of nearest neighbors with Mash/MinHash, clustering of homologous proteins with OrthoMCL, alignment of conserved clusters with MUSCLE, trimming with Gblocks, and concatenation followed by inference of a ML tree with RAxML. All complete Desulfovibrionaceae, Desulfobacteraceae, and Desulfobulbaceae genomes in the PATRIC database (n = 43, 21, and 11, respectively) were included in the analysis. The resulting tree was annotated using FigTree version 1.4.3. Subsequently, these three MAGs were compared to all publicly available Desulfovibrionaceae, Desulfobacteraceae, and Desulfobulbaceae genomes in GenBank (n = 122, 107, and 65, respectively) by average nucleotide identity (ANI) using fastANI, using >95% ANI as the intraspecies threshold and <83% as an interspecies threshold. Additionally, PATRIC’s “Protein Family Sorter” tool was used to identify sulfate-reducing proteins within the three SRM MAGs. To identify any polymer degradation protein sequences within the MAGs, pCDS were generated using Prodigal version 2.6.3 and were aligned against the custom hydrolase database (described above) using the same search parameters.

18Predicted coding sequences (pCDS) from each of the 12 metagenomes were aligned against three databases representing the three steps of the dissimilatory sulfate reduction pathway. The database for the first step of the pathway (reduction of sulfate to adenylyl sulfate [APS]) included all sulfate adenylyltransferase (SAT/MET3) protein sequences in the UniProtKB protein database. The database for the second step of the pathway (reduction of APS to sulfite) included all APS reductase (AprBA) protein sequences in the UniProtKB protein database. To address redundancy within both databases, CD-HIT version 4.7 was used to cluster protein sequences with default settings. This resulted in 612 and 239 protein sequences for the SAT/MET3 and AprBA databases, respectively. These two custom-made databases are also publicly available (see text footnote 1). For the third step of the pathway (reduction of sulfite to sulfide), metagenomes were aligned against a published database of DsrAB protein sequences. Alignments were performed using BLAST+ version 2.6.0 (blastp) with an e cutoff of 10−5, a minimum alignment length of 30 amino acids, and a minimum percent identity of 50%. Positive alignments were normalized to the total pCDS, and the abundance of SAT/MET3, AprBA, and DsrAB sequences was compared between all samples (n = 11; PET2 failed normality testing).

19Unless stated otherwise, R version 3.4.0 was used for statistical analysis of data. The combination of Shapiro–Wilk tests and quantile-quantile plots were used to test data for normality. One-way ANOVAs were generated using the R package multcomp, using a Tukey’s post hoc test with Westfall values.

20All sequence reads were made available through the project accession ERP017130 at the European Nucleotide Archive. The six MAGs were deposited as Madre1, 2, 3, 4, 5, and 6 at GenBank under the accession numbers QZKZ00000000, QZLA00000000, QZLB00000000, QZLC0000000, QZLD00000000, and QZLE00000000, respectively.

21Overall trends in the seawater-, ceramic-, PET-, and PHA-associated community compositions were visualized using the taxonomic rank of order based on the normalized proportion of OTUs assigned to each of the 12 samples (three seawater, three ceramic, three PET, and three PHA samples; Figure 1). The comparison of weighted UniFrac values revealed a significant difference between all three biofilm communities (ceramic, PET, and PHA) and the seawater community (PERMANOVA; p <0.05). Additionally, a significant difference was found between the PHA biofilm communities and both the ceramic and PET biofilm communities (PERMANOVA; p <0.05). By contrast, there was no significant difference between the PET and ceramic biofilm communities. Hierarchical clustering based on a Bray-Curtis resemblance matrix of OTU taxonomic assignments determined that the PET and ceramic biofilm communities were more closely related to the PHA biofilm community than to the seawater community (PERMANOVA; p <0.05).

22At the phylum level, OTUs belonging to Proteobacteria were the most abundant among all four community types, representing 42, 55, 59, and 72% of all OTUs in seawater, ceramic, PET, and PHA samples, respectively (Supplementary Table S2). Cyanobacteria was the second most abundant phylum in the seawater community (22%), making up a larger proportion of the community compared to the PET and ceramic biofilm communities (~6% for each), and even more so compared to the PHA biofilm community (0.7%). The presence of Chloroflexi (4%), Spirochaetes (4%), and Firmicutes (2%) among the most abundant phyla was unique to the PHA biofilm communities. A detailed report of taxonomic rank based on normalized proportion is provided in Supplementary Table S2.

23The OTUs assigned to orders Subsection I of Cyanobacteria, Oceanospirillales, SAR11, and Sphingobacteriales were among the top five most abundant in the seawater community only, making up 20, 11, 8, and 7% of the total OTUs in seawater, respectively. The orders Desulfobacterales, Rhodobacterales, Xanthomondales, and Cellvibrionales were all among the top five most abundant orders in both the PET and ceramic biofilm communities. The abundance of Proteobacteria (72% of all OTUs) within the PHA biofilm community reflected the dominance of the order Desulfobacterales, an order of SRM that represented 37% of all OTUs and more than half of all Proteobacteria. Similar to the trends with phyla, the abundant orders within the PHA biofilm community were unique, with the orders of Gammaproteobacteria incertae sedis, Anaerolineales, and Spirochaetales being in the top five most abundantly assigned OTUs in the PHA biofilm communities.

24Two genera of Cyanobacteria (Synechococcus and Prochlorococcus) represented almost 20% of all OTUs within the seawater community but made up a very small proportion (<1%) of all three biofilm communities. Members of uncultured genera from Desulfobacteraceae, Rhodobacteraceae, and Flammeovirgaceae were among the five most abundant genera in the PET and ceramic biofilm communities, representing almost 20% of all OTUs. Three genera within Desulfobacteraceae and an uncultured genus of Desulfobulbaceae combined to represent 20% of all OTUs within the PHA biofilm, reinforcing the dominance of SRM in that community (Supplementary Table S2).

25The comparison of normalized sequence abundance data revealed a significant enrichment of depolymerase and esterase gene sequences within the PHA biofilm communities versus the seawater, ceramic, and PET biofilm communities (ANOVA; p <0.05; Figure 2). The increase in depolymerases was nearly 20-fold, while the increase in esterases was approximately 2.5-fold. The significant increase in depolymerases was largely attributed to a polyhydroxybutyrate (PHB) depolymerase that constituted 60% of all hydrolases in the PHA biofilm communities but less than 0.4% of all hydrolases in other community types (Figure 3). Conversely, no hydrolases were enriched in the PET biofilm community versus either the ceramic biofilm or seawater communities. Hierarchical clustering based on a Bray-Curtis resemblance matrix of hydrolase relative abundance confirmed that PHA biofilm enzyme pools were significantly different from all other community types (SIMPROF; p <0.05). Seawater enzyme pools were also unique, while the PET and ceramic biofilms were indistinguishable (Figure 3).

26To characterize the enriched PHB depolymerases within the PHA biofilm communities, the three PHA metagenomes were co-assembled, resulting in 118,520 contigs totaling 245 Mbp. The N50 of the co-assembly was 2,119 bp, and the lengths of the smallest and longest contigs were 1,000 and 107,130 bp, respectively. The average contig length was 2,251 bp. The PHB depolymerase gene sequence was detected in 46 contigs. The average sequence identity between the 46 PHB depolymerase gene sequences was 54.0% with a range of 7.8–98.9%, and the ML tree illustrated the diversity of those 46 genes along with PHB depolymerase gene sequences from 10 known PHB degraders (Supplementary Figure S3).

27Six high-quality MAGs were recovered from the co-assembled PHA biofilm metagenomes. Initial MAG identification was inferred with CheckM: Madre1 was placed in the genus Desulfovibrio, Madre2 in the Desulfobacteraceae family, Madre3 in the Desulfobulbaceae family, Madre4 in the Spirochaetaceae family, and Madre5 and Madre6 within the Gammaproteobacteria order (Table 1). The three SRM MAGs (Madre1, 2, and 3) averaged 391 contigs, a N50 value of 18,335 bp, and a maximum contig length of 76,914 bp. The genome sizes ranged from 3,852,664 to 5,613,049 bp (Madre1 and Madre2, respectively), while the number of genes ranged from 3,829 to 5,179 (Madre1 and Madre2, respectively), and the GC content ranged from 46.9 to 59.8% (Madre3 and Madre1, respectively).

28The identities of the three PHA-associated SRM genomes were further explored by building a genome-scale phylogenetic tree that included 43 Desulfovibrionaceae, 21 Desulfobacteraceae, and 11 Desulfobulbaceae reference genomes (Figure 4). The phylogeny confirmed the initial CheckM-based placement of the three MAGs and further refined their closest relatives. Madre1 was related to the type strain Desulfovibrio gigas DSM 1382 (ATCC 19364). Madre2 was related to the uncultured Desulfobacula toluolica Tol2, an aromatic carbon-degrading SRM isolated from a seawater pond in Massachusetts. Madre3 was related to the uncultured Desulfofustis sp. PB-SRB1 that was recovered from the “pink berry” consortia of the Sippewissett salt marsh. However, a species could not be assigned to Madre 1, 2, or 3 given that each showed less than 83% ANI with the 122 Desulfovibrionaceae, 107 Desulfobacteraceae, and 65 Desulfobulbaceae reference genomes available in Genbank. The full tree, without collapsed branches, is included in the supplemental material (Supplementary Figure S4).

29Analysis of the three SRM MAGs revealed the presence of SAT/MET3 (Madre1, 2, and 3), AprBA (Madre2, and 3), and DsrAB (Madre1, 2, and 3) protein sequences. Additionally, the alignment of their pCDS to the custom hydrolase database revealed the presence of two distinct PHA depolymerase sequences in both Madre1 and Madre3.

30The enrichment of dissimilatory sulfur reduction protein sequences within the PHA biofilm communities was assessed through the alignment of metagenomic pCDS against databases of sulfate adenylyltransferases (SAT/MET3), adenylyl sulfate (APS) reductases (AprBA), and sulfite reductases (DsrAB). The comparison of normalized sequence abundance revealed no significant differences in SAT/MET3 protein abundances. A significant enrichment in AprBA and DsrAB was detected in the PHA biofilm communities (ANOVA; p <0.05; Figure 5).

31Plastic is a prominent marine pollutant. Numerous studies have characterized microbe-plastic interactions in pelagic systems, yet the effects of plastic pollution on benthic microbiota and biogeochemical cycling remain unclear. To address this knowledge gap, we conducted an in situ microcosm designed to ask how the benthic microbial community responds to plastic and bioplastic in coastal marine sediments. Coastal sediments are especially vulnerable to plastic loading due to their proximity to population centers and urban environments, and the majority of plastic debris entering the oceans accumulates in coastal zones.

32The comparison between surface-associated and free-living lifestyles has been a long-studied question in microbial ecology and, thereby, considerable literature is available to inform the study of surface-associated bacterial communities. Based on this literature, it was hypothesized that all surface-associated communities (i.e., ceramic, plastic, and bioplastic) would be distinct in comparison to free-living communities, and the results of this study supported that hypothesis. The more salient question was the distinction between different surfaces, and this study revealed that PET was not colonized by a distinct community in comparison to the ceramic biofilm control. This finding is supported by a recent study demonstrating that pelagic microbial communities associated with PET were indistinguishable from those associated with a glass biofilm control. Similarly, it was previously demonstrated that inert surfaces such as glass, ceramic, and coral skeleton have little influence on marine microbial community composition. The lack of distinction between PET and ceramic could also be due to the remote location of the study site seeing that previous reports of plastic degradation have occurred in highly polluted habitats that would select for and enrich plastic degraders.

33Whereas PET was not colonized by a distinct microbial community, the introduction of PHA promoted a significant and distinct response. In particular, SRM were the dominant members of the PHA-associated assemblage (see Figure 3). Typically, the three most common SRM families (Desulfobacteraceae, Desulfobulbaceae, and Desulfovibrionaceae) account for between 5 and 20% of the total bacterial community in estuarine sediments. In this study, these three families made up a similar fraction of the plastic- and ceramic-associated benthic communities (9 and 12%, respectively) but they accounted for over 39% of the PHA-associated benthic community.

34Naturally occurring PHAs are carbon and energy storage polymers that are accumulated as intracellular granules that aid survival during periods of nutrient limitation and environmental stress. The ability to produce and degrade PHA is thought to be widespread but its prevalence among SRM and its importance in carbon and sulfur cycling has not been thoroughly characterized. A previous study has reported PHA accumulation in Desulfococcus multivorans, Desulfobotulus sapovorans, Desulfonema magnum, and Desulfobacterium autotrophicum. Importantly, two culture-dependent studies have shown that addition of PHA to incubation bottles was correlated with sulfide production in anoxic lake sediments. Taken together, the accumulation of PHA in SRM, the PHA stimulation of sulfide production, and the PHA-selection of SRM (in this study) indicate that PHA is an important carbon source for sulfate reduction. However, it remains unknown if the sulfate reducers in our study can directly degrade PHA or if they rely on primary fermenters to first degrade the polymer into simple organic substrates.

35A key reaction in PHA degradation is its depolymerization. In this study, a nearly 20-fold increase in the abundance of depolymerases indicated that PHA stimulated the growth of PHA-degrading bacteria. Further, the large diversity of the 46 depolymerases, recovered from the metagenome-assembled genomes, suggested these enzymes were distributed across a mixed microbial consortium. Although we did not measure enzyme activity, the increased abundance and diversity of these depolymerases support the hypothesis that PHA biofilms were sites of enhanced enzyme activity. Additional support for this hypothesis will be made available in a parallel and forthcoming 15-month study wherein weight loss and scanning electron microscopy were used to quantify PHA degradation at the same study site (Figure 6).

36Previous studies have shown, using adenylyl sulfate reductases (aprBA) and dissimilatory sulfite reductases (dsrAB) as functional markers, that sulfate reduction is more diverse and widespread than previously thought. Further, the detection of apr and dsr genes in new phyla, made possible through the construction of genomes from metagenomes, has expanded the known diversity of sulfate-reducing microorganisms. Here, we demonstrated that bioplastic communities were significantly enriched with both sulfate and sulfite reductases. We also reconstructed genomes from metagenomes to discover new Desulfovibrio, Desulfobacula, and Desulfofustis species (Madre 1, 2, and 3, respectively). While it remains unknown if SRM possessing depolymerases (i.e., Madre1 and Madre 3) can directly degrade PHA, an ongoing study will shed light on their growth requirements. Despite this open question, the increased abundance of sulfate-reducing enzymes and the discovery of three uncultured sulfate-reducing species lend additional support to the hypothesis that sulfate reduction in marine sediments is stimulated by the addition of PHA.

37In a scenario where bioplastic use and pollution are more common, what remains unclear is the larger effect of the PHA-selection for SRM. Sulfate-reducing bacteria play a critical role in biogeochemical cycling as sulfate reduction is responsible for ~50% of organic carbon mineralization in marine sediments, and the stimulation of SRM has been shown to suppress the growth of methanogenic archaea through the diversion of carbon flow from methane to carbon dioxide. For example, the stimulation of SRM in coastal rice fields, through the application of sulfate-containing amendments, was shown to inhibit methane production. It is thus possible that bioplastic pollution and its stimulation of SRM could have unintended consequences that affect the balance between sulfate reduction and methanogenesis. More generally, it is also possible that bioplastic loading could alter the natural syntrophic cycling of PHA in marine sediments.

38In conclusion, this study demonstrated that the introduction of plastic had no measurable effect on the benthic microbial community. By contrast, the introduction of bioplastic selected for a distinct microbial community that was enriched for hydrolases and dominated by SRM. Recovered genomes, representing the three most common families of SRM, contained depolymerases as well as aprBA and dsrAB reductases. Findings indicate that SRM play an important role in PHA degradation in coastal marine sediments. Given that sulfate reduction is a key process in the oceanic sulfur cycle, we recommend that future scientific investigation, government legislation, and best-management practices related to plastic pollution consider the effects of plastic as well as bio-based alternatives on benthic biogeochemical activities.

39The datasets generated for this study can be found in ENA/EBI Metagenomics & NCBI, ERP017130 & QZKZ00000000, QZLA00000000, QZLB00000000, QZLC0000000, QZLD00000000, QZLE00000000.

40All authors contributed to this study’s design, implementation, analyses, and writing.

41The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.