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

Plastic is incredibly abundant in marine environments but little is known about its effects on benthic microbiota and biogeochemical cycling. This study reports the shotgun metagenomic sequencing of biofilms fouling plastic and bioplastic microcosms staged at the sediment–water interface of a coastal lagoon. Community composition analysis revealed that plastic biofilms were indistinguishable in comparison to a ceramic biofilm control. By contrast, bioplastic biofilms were distinct and dominated by sulfate-reducing microorganisms (SRM). Analysis of bioplastic gene pools revealed the enrichment of esterases, depolymerases, adenylyl sulfate reductases (aprBA), and dissimilatory sulfite reductases (dsrAB). The nearly 20-fold enrichment of a phylogenetically diverse polyhydroxybutyrate (PHB) depolymerase suggests this gene was distributed across a mixed microbial assemblage. The metagenomic reconstruction of genomes identified novel species of Desulfovibrio, Desulfobacteraceae, and Desulfobulbaceae among the abundant SRM, and these genomes contained genes integral to both bioplastic degradation and sulfate reduction. Findings indicate that bioplastic promoted a rapid and significant shift in benthic microbial diversity and gene pools, selecting for microbes that participate in bioplastic degradation and sulfate reduction. If plastic pollution is traded for bioplastic pollution and sedimentary inputs are large, the microbial response could unintentionally affect benthic biogeochemical activities through the stimulation of sulfate reducers.


INTRODUCTION
Microorganisms rapidly colonize and form biofilms on biotic and abiotic surfaces in marine environments (Dang and Lovell, 2016). As such, marine microbial communities are commonly distinguished as free-living or particle-associated (DeLong et al., 1993;Crump et al., 1999;Hollibaugh et al., 2000). A particle-associated lifestyle is thought to provide microorganisms with increased nutritional resources and environmental stability (Azam et al., 1994). Detrital aggregates, in particular, are hot spots of microbial diversity that are fundamentally distinct in comparison to free-living communities (DeLong et al., 1993;Crump et al., 1999). Particles are also habitats of enhanced enzyme activity (Karner and Herndl, 1992;Kellogg and Deming, 2014) that contribute to the rapid turnover of particulate organic matter (Biddanda and Pomeroy, 1988;Azam and Long, 2001) and it is clear that particle-associated microorganisms play an important role in biogeochemical processes (Ploug et al., 1999;Seymour et al., 2017).
Plastic pollution has introduced a new surface for microbial colonization and biofilm formation (Zettler et al., 2013;Oberbeckmann et al., 2014). 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 (Geyer et al., 2017). Oceans serve as a sink for said waste with an estimated 4.8-12.7 million metric tons entering the oceans annually (Jambeck et al., 2015;Geyer et al., 2017). 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 (Bryant et al., 2016;Oberbeckmann et al., 2016).
Several studies have characterized pelagic microbe-plastic interactions (Carson et al., 2013;Zettler et al., 2013;Bryant et al., 2016) but relatively few studies have characterized microbe-plastic interactions in benthic systems (Harrison et al., 2014;De Tender et al., 2015;Nauendorf et al., 2016). In the North Pacific Gyre, a previous study reported that microplastics were colonized by diatoms (Carson et al., 2013) and a separate study reported an increased concentration of Chl a and an increased abundance of nitrogenase genes in microplastic biofilms (Bryant et al., 2016). 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 (Barnes et al., 2009;Andrady, 2011). Sinking is therefore an important process and benthic systems are sinks for plastic debris accumulation (Thompson et al., 2004;Van Cauwenberghe et al., 2015).
In coastal regions, benthic microorganisms support primary producers in the overlying water column by the remineralization of organic-rich debris (Klump and Martens, 1981). The remineralization of carbon in anoxic coastal sediments is predominantly carried out via sulfate reduction (Jørgensen, 1977). 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 (Vigneron et al., 2018) and recent studies have clearly shown that sulfate and sulfite reducers are more abundant, diverse, and widespread than previously believed (Müller et al., 2014;Anantharaman et al., 2018;Hausmann et al., 2018;Vigneron et al., 2018;Zecchin et al., 2018). Yet, the effect plastic pollution may have on sulfate-reducing microorganisms (SRM) is an open question.
In 2018, the United Nations Environmental Program reported that more than 60 countries approved bans or levies on single-use plastics (UNEP, 2018). 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 (European Commission, 2018). 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 (Aeschelmann and Carus, 2015). 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 (Eich et al., 2015;Pauli et al., 2017), little is known about how microorganisms respond to bioplastic versus petroleum-based plastic.
The 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.

Site Description
The 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 (Tunnell, 2002). The dredge material island itself is remote to the nearest population center (Corpus Christi, TX) and is therefore relatively free of plastic pollution.

Sample Deployment and Collection
A 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 (Andrady, 2011). 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 ).

DNA Isolation
Genomic 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 (Zhou et al., 1996). 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 (A 260 /A 280 and A 260 /A 230 ) 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.

Metagenome Sequencing
Metagenomic 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 (Mitchell et al., 2016) version 3.0, which includes an automated workflow for read processing. Briefly, overlapping reads were merged with SeqPrep (John, 2011), low-quality reads and adapter sequences were trimmed using trimmomatic (Bolger et al., 2014), and reads less than 100 bp were removed using BioPython (Cock et al., 2009). Identification of the 16S rRNA reads was performed using rRNASelector (Lee et al., 2011) and FragGeneScan (Rho et al., 2010) 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.

Microbial Community Composition
Operational taxonomic units (OTUs) were assigned using the QIIME (Caporaso et al., 2010) version 1.9 closed-reference OTU picking protocol and the SILVA 128 SSU 97% database (Quast et al., 2013) with reverse strand matching enabled. Beta-diversity was analyzed using weighted UniFrac (Lozupone et al., 2011) 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 (Clarke and Gorley, 2015) (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.

Enrichment of Alpha/Beta Hydrolases
Predicted coding sequences (pCDS) from each of the 12 metagenomes were aligned against a modified database 1 of protein sequences from the ESTHER database of alpha/beta hydrolases (Hotelier et al., 2004). Proteins included in the modified database were restricted to depolymerases, esterases, lipases, and cutinases as these were previously implicated in plastic-microbe interactions (Yoshida et al., 2016). To reduce redundancy within the database, CD-HIT (Clarke and Gorley, 2015) 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+ (Camacho et al., 2009) 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) (Clarke et al., 2008) 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).

Relatedness of Polyhydroxybutyrate Depolymerases
Processed reads for each of the PHA samples (n = 3) were co-assembled using MEGAHIT (Li et al., 2015) version 1.1.1. Co-assembly was performed using the "meta-large" preset for large and complex communities, and MEGAHIT was called as follows: megahit -presets meta-large -r input --min-contig-len 1,000 -o output -t 40.
To 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 (Nguyen et al., 2015) version 1.6.1. Sequences were aligned using the M-Coffee mode of T-Coffee (Notredame et al., 2000) that combines results from eight different aligners and the tree was generated with 1,000 ultrafast bootstraps (Minh et al., 2013) using the best-fit model (TPM3u + F + R3) as determined by ModelFinder (Kalyaanamoorthy et al., 2017). The phylogenetic tree was annotated with FigTree version 1.4.3 2 . Sequence similarity between the 46 PHB depolymerase gene sequences was determined using trimAl (Capella-Gutiérrez et al., 2009)

Recovery and Analysis of Genomes of Sulfate-Reducing Microorganisms
To characterize individual SRM within the PHA biofilm community, metagenome-assembled genomes (MAGs) were recovered using a previously described method (Parks et al., 2017). Briefly, the processed reads were mapped to the co-assembled PHA metagenomic contigs (produced using MEGAHIT above) with BWA (Li and Durbin, 2010) version 0.7.15-r1142 using the BWA-MEM algorithm with default parameters. Genomes were recovered with MetaBAT (Kang et al., 2015) 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 (Parks et al., 2015) version 1.0.11 and RefineM version 0.0.23, as described previously (Parks et al., 2017). Of the 46 bins produced, six MAGs were retained for further exploration based on criteria 2 http://tree.bio.ed.ac.uk/software/figtree/ adopted by Parks et al. (Parks et al., 2017): 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.
The three SRM MAGs were analyzed with the Pathosystems Resource Integration Center's (PATRIC) (Wattam et al., 2017) comprehensive genome analysis service. The PATRIC database contains over 190,000 bacterial genomes and has been increasingly used in environmental studies (Harke et al., 2015;Skennerton et al., 2016;Ortiz-Álvarez et al., 2018).The automated PATRIC service includes annotation with RASTtk (Brettin et al., 2015), prediction of nearest neighbors with Mash/MinHash (Ondov et al., 2016), clustering of homologous proteins with OrthoMCL (Enright et al., 2002), alignment of conserved clusters with MUSCLE (Edgar, 2004), trimming with Gblocks (Talavera and Castresana, 2007), and concatenation followed by inference of a ML tree with RAxML (Stamatakis, 2014). 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 (Jain et al., 2018), 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 (Hyatt et al., 2010) version 2.6.3 and were aligned against the custom hydrolase database (described above) using the same search parameters.

Dissimilatory Sulfur Reduction Potential
Predicted 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 (Clarke and Gorley, 2015) 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 (Müller et al., 2014) protein sequences. Alignments were performed using BLAST+ (Camacho et al., 2009) 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).

Statistical Analyses
Unless stated otherwise, R (R Core Team, 2017) 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.

Data Availability
All 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.

Microbial Community Composition
Overall 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).
At 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.
The 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.
Two 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).

Enrichment of Alpha/Beta Hydrolases
The 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).

Relatedness of Polyhydroxybutyrate Depolymerases
To 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).

Recovery and Analysis of SRM Genomes
Six 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 FIGURE 2 | Bar plot showing the relative abundance of depolymerases, esterases, lipases, and cutinases. Gene sequences were normalized to the total number of pCDS in each sample, with individual samples plotted and grouped together based on community association (seawater, ceramic, PET, and PHA). Significance was determined using ANOVA (***p = <0.001, n = 3).
The 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) (Le Gall, 1963;Morais-Silva et al., 2014). Madre2 was related to the uncultured Desulfobacula toluolica Tol2, an aromatic carbon-degrading SRM isolated from a seawater pond in Massachusetts (Wöhlbrand et al., 2013). Madre3 was related to the uncultured Desulfofustis sp. PB-SRB1 that was recovered from the "pink berry" consortia of the Sippewissett salt marsh (Wilbanks et al., 2014). 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).
Analysis 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.

Dissimilatory Sulfur Reduction Potential
The 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).

DISCUSSION
Plastic 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 (Barnes et al., 2009), and the majority of plastic debris entering the oceans accumulates in coastal zones (Barnes et al., 2009;Thompson et al., 2009). The 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 (Biddanda and Pomeroy, 1988;Karner and Herndl, 1992;DeLong et al., 1993;Crump et al., 1999;Hollibaugh et al., 2000;Kellogg and Deming, 2014). 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 (Oberbeckmann et al., 2016). Similarly, it was previously demonstrated that inert surfaces such as glass, ceramic, and coral skeleton have little influence on marine microbial community composition (Witt et al., 2011). 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 (Nanda et al., 2010;Singh and Gupta, 2014;Yoshida et al., 2016).
Whereas 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 (Bowen et al., 2012;Cheung et al., 2018). 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.
Naturally 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 (Dawes and Senior, 1973;Obruca et al., 2018). The ability to produce and degrade PHA is thought to be widespread (Reddy et al., 2003) 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 (Hai et al., 2004). Importantly, two culturedependent studies have shown that addition of PHA to incubation bottles was correlated with sulfide production in anoxic lake sediments (Mas-Castellà et al., 1995;Urmeneta et al., 1995). 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.
A 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).
Previous 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 (Müller et al., 2014;Anantharaman et al., 2018;Hausmann et al., 2018). 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 Hausmann et al., 2018). 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.
In 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 (Jørgensen, 1982), and the stimulation of SRM has been shown Frontiers in Microbiology | www.frontiersin.org to suppress the growth of methanogenic archaea through the diversion of carbon flow from methane to carbon dioxide (Gauci et al., 2004). For example, the stimulation of SRM in coastal rice fields, through the application of sulfate-containing amendments, was shown to inhibit methane production (Lindau et al., 1994;Denier van der Gon et al., 2001). 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.
In 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 bestmanagement practices related to plastic pollution consider the effects of plastic as well as bio-based alternatives on benthic biogeochemical activities.

AUTHOR CONTRIBUTIONS
All authors contributed to this study's design, implementation, analyses, and writing.

ACKNOWLEDGMENTS
This study would not have been possible without the support of our friends and colleagues. Specifically, we thank Thomas Merrick for his assistance with data analysis, Robert "Bobby" Duke and the entire Center for Coastal Studies staff for assistance in the field, and Nicole C. Elledge for her assistance with large-scale ANI comparisons.
TABLE S2 | The top five members from each taxonomic rank based on normalized proportion per sample. Means (n = 3) plus or minus the standard error of mean are reported for each of the taxonomic assignments. FIGURE 6 | Scanning electron microscopy image of PHA pellets used to visualize signs of biodegradation. Panel (A) displays a pellet that was not exposed to benthic microbial communities, while Panel (B) shows a pellet exposed for 230 days.
deployed at the water-sediment interface using a randomized experimental design as shown in panel (D).