Introduction

An estimated total of 7 billion metric tons of plastic waste has now been produced globally while approximately only 9% of this is recycled and 79% has been discarded in landfills or the environment [1]. In marine and other aquatic environments, plastics cause a range of negative environmental impacts: directly, through entanglement [2, 3] or ingestion [4,5,6,7,8], or indirectly, through the transfer of toxic chemicals [9, 10]. Despite that the majority of marine plastic waste originates from land, less is known about the impacts of plastic waste in terrestrial ecosystems [11], and there are numerous factors that can affect the fate, transport, and impacts of plastics in all biospheres.

These factors include transport by air [12], rain [13], rivers [14] and currents [15], (de)sorption of chemicals [16], photo- or mechanical degradation or fragmentation [17] and microbial colonization and possibly degradation [18]. Deleterious impacts of plastics [6], their potential to transport invasive species across entire ocean basins [3, 19,20,21] and their biodegradation by isolated bacteria and fungi have been studied for over half a century [22]. The microbial communities colonizing plastics—commonly termed the “Plastisphere” [23]—however, have only been specifically investigated more recently. A call for research into the interactions between microorganisms and plastics at the beginning of 2011 marks almost a decade of Plastisphere research [24]. On 5 January, 2020 there were 50 publications that studied the Plastisphere using Illumina Next Generation Sequencing (NGS) methods. Although this is an emerging area of research, most Plastisphere studies have been focussed on the marine environment.

Microbial members of the Plastisphere have been found to be: (i) different from communities that colonize other surfaces [25,26,27]; (ii) not different from communities that colonize other surfaces [28]; (iii) only different from communities colonizing other surfaces under specific environmental conditions [29, 30] or at specific time points [31]; (iv) more diverse than other microbial communities [32]; (v) less diverse than other microbial communities [33, 34]; (vi) potentially degrading the plastics that they colonize [35, 36]; (vii) capable of degrading plastic additives [33, 37, 38]; and (viii) pathogenic and/or carrying antimicrobial resistance genes [39,40,41,42]. The marine Plastisphere has been heavily reviewed within the last year, e.g., [18, 43,44,45,46], and there are also several recent reviews on plastic biodegradation, e.g., [47,48,49]. However, definitive answers on the metabolic capabilities of the Plastisphere or the factors that drive its formation and composition are unknown. Indeed, a re-analysis of a small subset of Plastisphere studies (n = 5) by Oberbeckmann and Labrenz [45] revealed that salinity along with other environmental factors appeared to have a larger effect on community composition than substrate.

To date, no large-scale meta-analysis of the Plastisphere has been conducted, despite that all NGS data collected are theoretically comparable. This is likely because these data are not always deposited in publicly accessible databases (as best practice would dictate). In addition, the use of different methods for processing these data means that direct comparisons have not been possible without substantial bioinformatic analyses. In this study, we conducted a re-analysis of all studies from 2010 to 2019 (for which sequencing data were already accessible or were made available upon request) that use Illumina NGS to characterize the Plastisphere using the 16S rRNA gene. We aimed to determine whether taxa identified as potential plastic biodegraders or potential pathogens -and that were higher in abundance on plastics than other samples- were significant across multiple studies and environments. We then investigated which environmental and methodological factors were shaping the Plastisphere. We classified all sequences from these studies to amplicon sequence variants (ASVs) and used a phylogeny-based method [50] to overcome the problems presented by the use of different primer pairs. This allowed us to identify the common taxa between these studies and to use random forest models to draw conclusions on the over-arching factors that shape the Plastisphere.

Materials and methods

Experimental design

A literature search was performed on 5 January, 2020 using the search terms “Plastics plastisphere”, “Plastics microbial community”, and “Plastics microbial degradation” in both the Web of Science Core Collection and Science Direct (Supplementary Table 1A). The search was limited to studies that fit the following criteria: (i) were published between 2010 and 2019; (ii) had original data (iii) characterize the biofilm formed on nonbiodegradable plastics; and (iv) use Illumina NGS. This resulted in 50 studies for inclusion in this meta-analysis (Supplementary Table S1B); 41 studies that characterized only the 16S rRNA gene [25, 27,28,29,30,31,32,33, 35, 38, 41, 42, 51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79], two studies that characterized only the 18S rRNA gene [34, 80], two studies that characterized the 16S and 18S rRNA genes [26, 81], two studies that characterized the 16S and ITS2 rRNA genes [82, 83], two studies that used shotgun metagenomics [84, 85], and one study that used shotgun metagenomics and characterized the 16S rRNA gene [86].

Of these studies, 34 had sequencing data that were already publicly accessible. Requests for raw sequencing data and metadata were made to the corresponding and first authors of the remaining 16 studies, resulting in the provision of these data for a further five studies. This resulted in 39 studies with datasets that were available for inclusion: 31 that characterized only the 16S rRNA gene [25, 27, 28, 30,31,32,33, 35, 38, 52, 53, 55,56,57,58,59, 61, 62, 64,65,66, 68, 69, 71,72,73,74,75,76, 78, 79], two that characterized only the 18S rRNA gene [34, 80], two that characterized the 16S and ITS2 rRNA genes [82, 83], two that used shotgun metagenomics [84, 85], one that characterized the 16S and 18S rRNA genes [26], and one that used shotgun metagenomics and characterized the 16S rRNA gene [86] (Supplementary Table S1B). Relatively few studies focussed on the 18S and ITS2 rRNA genes and these used primers that targeted different regions and were from different environments. We therefore focussed on the 16S rRNA gene here, meaning that we included a total of 35 studies [25, 27, 28, 30,31,32,33, 35, 38, 52, 53, 55,56,57,58,59, 61, 62, 64,65,66, 68, 69, 71,72,73,74,75,76, 78, 79, 82, 83, 86] with 2,229 samples between then. Publicly accessible data containing sequencing reads for the remaining studies were downloaded (primarily from the NCBI SRA database) and, if necessary, files were converted to match a format that was compatible with QIIME2 (i.e., GNU zipped FASTQ files). Additional requests for raw sequencing data were made to the authors of studies where the forward and reverse reads were already joined, or the primers were already removed (full details can be found in Supplementary Table S1B). All available metadata were collected and supplemented with any additional information present in the supplementary information of published studies (Supplementary Table S2). Where details of salinity and temperature were not given, these were estimated based on typical characteristics in these areas at the time of year the samples were collected. If no light regime was specified, ambient light was assumed. If sample names were not given that could be matched between supporting information of the paper and metadata given to the NCBI SRA or plastic type was not determined, these were classified as “unknown plastic”.

We added metadata categories that were not in any of the original studies: (i) environment—each study was classified as terrestrial, marine, freshwater, or aquatic (e.g., the environment was brackish or the experiment was carried out in a different system, such as in mariculture); (ii) whether the study was carried out in a laboratory or in the field; (iii) incubation or collection—whether the plastics used were incubated for a known length of time or were collected after an unknown residence time; (iv) general incubation time—samples were classified as early (≤7 days incubation), late (>7 days incubation), or collection (samples were collected after an unknown residence time); (v) water or sediment—the plastics were from/were incubated either in the water column or on the sea floor/in soil in the terrestrial environment; (vi) source—the material used for microbial community characterization was classified to the general categories of plastic, water, sediment, organic, not plastic (i.e., an inert control surface, such as glass or metal), other or blanks and positive sequencing or methodological controls; and (vii) general plastic type—samples were classified as aliphatic (i.e., PE or PP), other plastic (i.e., plastics that contain other functional groups, e.g., PET, PS, and PVC), unknown plastic (the plastic type was not determined), biofilm (the sample was from a control substrate, such as glass or a leaf), planktonic, or blank (i.e., sequencing or methodological controls).

16S amplicon sequence processing (per study)

Data processing followed the standard operating procedure suggested in Comeau et al. [87] and used QIIME2 (2019.10 core distribution [88]). Raw forward and reverse read files were imported to the QIIME2 format after an initial visualization of read quality using the packages FastQC (Babraham Bioinformatics) and MultiQC [89]. Cutadapt [90] was used to remove primers from reads and VSEARCH [91] was used to join paired-end reads. These steps were omitted for samples where primers were already removed or reads were already joined, respectively, and for one study that used three different reverse primers [59]. All low-quality reads were then filtered using default quality thresholds before using Deblur [92] to denoise sequences and resolve ASVs. Deblur was run using trim lengths determined by the read quality for each study. For Frére et al. [59], the ends of reads corresponding to the lengths of primer sequences were also trimmed. The forward reads only were used for several studies where low reverse read quality led to too few reads remaining after running Deblur and for which the forward reads were of high enough quality to be used (full details of the processing steps carried out for each study can be found in Supplementary Table S2).

Combined processing

All studies were merged using QIIME2’s merge and merge-seqs commands, then classified taxonomically using a classifier trained on the full-length 16S rRNA gene SILVA v132 database [93]. Classified sequences were filtered to remove mitochondria, chloroplasts, those that were unclassified at the kingdom level and those present at a cumulative abundance of ten or fewer. This left a median of 20,284 reads per sample (minimum 2 and maximum 995,391). Samples with <2,000 reads were removed, leaving 2,056 samples and 34 studies, and phylogenetic trees were built using SEPP with a reference phylogeny created using the SILVA v128 database [50].

Custom scripts that wrapped all commands were used to carry out further analyses in R (version 3.6.1) and Python (version 3.8.3) using the data exported from QIIME2 and the packages Biopython [94], csv, itertools, lifelines [95], math, matplotlib [96], numpy, os, pandas [97], pdf2image, pickle, scipy [98], scikit-bio, scikit-learn [99], and sinfo for Python as well as ape [100], compositions [101], dplyr [102], exactRankTests [103], ggnewscale [104], ggplot2 [105], ggtree [106], knitr [107], metacoder [108], microbiome [109], nlme [110], philr [111], phyloseq [112], reticulate [113], and vegan [114] for R. Several normalization methods were used to address the large disparity in the sequencing depths of different samples: (i) samples were rarefied to 2,000 and converted to relative abundance; (ii) samples were converted to relative abundance; (iii) samples were converted to a log scale (with a pseudo count of one); or (iv) samples were converted to a centered log ratio (CLR; with a pseudo count of half of the minimum nonzero count). ASVs with a maximum number of reads below 1% of the median number of reads per sample were removed—below 20 reads for (i) and 202 reads for (ii), (iii), and (iv)—and sequences were agglomerated at a height of 0.1 based on the SEPP insertion tree. This resulted in 4,469 ASVs remaining for (i) and 12,635 ASVs remaining for (ii), (iii), and (iv). For (i), (ii), and (iii), weighted and unweighted uniFrac [115] distances were calculated between all samples while for (iv), the log-ratio transformed data, Aitchison distances (i.e., Euclidean distances of CLR-transformed data) [116] as well as Phylogenetic Isometric Log-Ratio Transformation [111] distances were calculated between all samples. Unless otherwise mentioned, all analyses used the rarefied data.

Random forest model construction

To determine which taxa were most associated with groupings within each metadata category (e.g., environment, plastic type, salinity, incubation time; n = 20 categories), feature selection was performed using random forest models (either classification or regression models for discrete or continuous data, respectively) built using scikit-learn, each with 10,000 estimators and 80 and 20% of samples being used for training and testing, respectively. All taxa were scaled to the maximum value of taxon abundance prior to building the models. These were carried out separately for the taxonomic levels of phylum, class, order, family, genera, species, and ASV. Where groupings within metadata categories for samples were not known, these samples were removed from the model construction. To investigate the effect of normalization method on the selected features, these models were constructed for each of the four normalization methods described above. To determine the taxa that were most associated with different plastic types, random forest classification models (built using the same parameters as above) were constructed separately for each environment (marine, aquatic, freshwater or terrestrial) for samples grouped to general plastic type (as described above). Samples were not further grouped into laboratory and field studies because not all environments included both. These were again constructed for the taxonomic levels of phylum, class, order, family, genus, species, and ASV and for each of the four normalization methods described above. The classification accuracy for each random forest model is defined as the percentage of the time that the model can classify a sample to the correct grouping (e.g., as marine within the environment metadata category) and the feature importance is defined as the proportion that the classification accuracy would decrease without that taxon (feature) present. We calculated concordance (using the Python lifelines package) [95] between the classification accuracies as well as between mean feature importance values to assess the similarity of the results obtained by the random forest models constructed using different normalization methods.

Tests for differential abundance of taxa

ASVs shared between treatments were calculated for samples grouped by the environment. Differences between samples at early (≤7 days) or late (>7 days) incubation times were determined using Wilcoxon Rank Sum tests within the metacoder R package with holm-bonferroni false discovery rate correction. ANCOM tests for differential abundance (with holm-bonferroni false discovery rate correction) were performed on all taxa (separately for each taxonomic level) for the taxa normalized using rarefaction.

Results

Summary of included studies and sequences

In this study, we reanalyzed the 16S rRNA gene amplicon sequencing data obtained from 35 studies (Supplementary Section 1). All Plastisphere studies to date have examined data collected in the Northern hemisphere and most of the included studies were conducted in and around Europe (Fig. 1). After removing all samples with below 2,000 reads there were 2,056 samples remaining; 1,185, 316, 506, and 49 samples in the marine, freshwater, aquatic, and terrestrial environments, respectively. One study [64] was removed from further analyses because all samples (n = 9) had below 2,000 reads. The abundance of different substratum types depended on the environment and most samples were collected from the field after unknown environmental residence times (Fig. 1).

Fig. 1: Overview of the studies and samples included in the meta-analysis.
figure 1

Cumulative number of studies per year (A), study location (B), number of samples (C), relative abundance of sample type (D), and sample incubation time (E) for studies carried out in the marine, freshwater, other aquatic, and terrestrial environments. A, B Show all 50 studies whereas, C, D, and E show only studies/samples that were included. Studies for which data were not provided (neither publicly available nor provided upon request) are shown with transparent colors in A and white marker edges in B. Note that those studies shown for 2020 were already in press and available online by 5 January, 2020. See Supplementary Tables S1 and S2 for full details of all studies and samples included.

Full size image

Normalization method affects the biological interpretation of results

There was a large disparity in the number of reads and taxonomic richness per sample between different studies (Fig. S1). We therefore investigated several different methods for data normalization (further details are given in the methods section) on random forest classification accuracy and the concordance between taxa identified as important by each of these methods (Supplementary Section 2). These random forest models were constructed for: (i) all 20 metadata categories—including factors such as environment, geographic location, temperature, primer pair, and plastic type to varying degrees of specificity—for each taxonomic level: phylum, class, order, family, genus, species, and ASV (i.e., 560 random forest models were constructed in total, 140 for each normalization method); and (ii) general plastic type within each of the four environments (marine, aquatic, freshwater, and terrestrial) for each taxonomic level: phylum, class, order, family, genus, species, and ASV (i.e., 112 random forest models were constructed in total, 28 for each normalization method). For both sets of random forest models, the models constructed using rarefied data were on average the least accurate (61–67% and 53–67% mean classification accuracy across all metadata categories and all environments for general plastic type, respectively, for the taxonomic levels of phylum-ASV; Figs. S2 and S3) and the models constructed on data transformed to relative abundance were on average the most accurate (82–90% and 72–86% mean classification accuracy across all metadata categories and all environments for general plastic type, respectively, for the taxonomic levels of phylum-ASV; Figs. S2 and S3). Models constructed using the compositionally aware centered log-ratio transformed data were on average 1.1 and 10% (for all metadata categories and all environments for general plastic type, respectively) less accurate than models constructed on data transformed to relative abundance. Feature importance values are the proportion that the classification accuracy decreases without that feature. The features that have importance values of either above 0.01 or 0.005 are the same between all normalization methods across all taxonomic levels. Concordance in feature importance values (whether the ranking of features by their values is the same) was on average 0.94 across all taxonomic levels between the relative abundance, log and CLR-transformed data and 0.78 for the relative abundance, log and CLR-transformed data against the rarefied data (Supplementary Section 2).

Rarefying has previously been found to most effectively account for large differences in library sizes, including lowering the false discovery rate when there are large differences in library sizes [117]. We therefore compromised on lower random forest classification accuracy and used the rarefied data for the remainder of analyses. Most analyses presented here have also been conducted on the data normalized using the other methods, but these results are presented in Supplementary Section 2.

Diversity within and between studies

Beta diversity analyses (weighted and unweighted uniFrac distance) showed that there were significant differences between samples when they were grouped by study or environment (PERMANOVA and ANOSIM p = 0.001; Fig. 2). The largest differences between studies and environments were due to organisms that were present at lower abundances, as evidenced by higher pseudo-F and R statistics with the unweighted (pseudo-F = 100.42 and 55.149 and R = 0.247 and 0.897 for environment and study, respectively) than weighted (pseudo-F = 93.87 and 53.045 and R = 0.224 and 0.716 for environment and study, respectively) uniFrac distances (Figs. 2 and 3). Studies that were clearly less similar to other studies could be explained by them being laboratory-based [32, 71], focussed on anaerobic rather than aerobic communities [71], collecting samples from the deep sea [76] or sequencing amplicons that did not include the V4 16S rRNA gene region [52, 71]. Those that were particularly similar within [26] or between [72,73,74] studies could presumably be explained by very long incubation times (above 1 year), leading to community convergence, or having similar experimental setups and inoculums, respectively (Fig. 3).

Fig. 2: nMDS plots showing uniFrac distance between samples.
figure 2

nMDS plots showing weighted (A, B) or unweighted (C, D) uniFrac distance (i.e., accounting for taxon phylogeny with or without taxon abundance, respectively) calculated between all samples and shown in A and C as samples colored by environment and B and D as samples colored by study. Results of PERMANOVA and ANOSIM tests for significance between groups are shown on each plot.

Full size image
Fig. 3: Average similarity between samples within a study versus between studies.
figure 3

Average similarity (determined by weighted or unweighted uniFrac; i.e., accounting for taxon phylogeny with or without taxon abundance; top or bottom, respectively) between samples within a study versus between studies, with white boxes showing samples grouped by environment. Study names are colored by environment, with green, purple, blue, and orange being for marine, aquatic, freshwater, and terrestrial, respectively (as in Figs. 1 and 2).

Full size image

When samples are split to environment and substratum type (i.e., the groupings shown in Fig. 1D), they have a tendency to group by environment, although this is not true for all cases (Fig. 4). The Proteobacteria dominate (above 50% relative abundance) in all but the terrestrial environment, freshwater planktonic, and aquatic control biofilm groups, which have larger proportions of Actinobacteria and Planctomycetes, respectively. Medians of Simpson’s diversity indices were above 0.8 for all sample groups apart from the aquatic blank group (0.3), which had a high relative abundance of Vibrionales (i.e., above 15%) and was the only group with >1% of the extremophilic Deinococcus–Thermus phylum (i.e., almost 7%). The number of ASVs in each environment (only those that were >1% relative abundance are shown) was related to the number of samples from that environment, with there being the highest and lowest numbers of both ASVs and samples in the marine and terrestrial environments, respectively (Figs. 1 and 4), and each substratum type within each environment had ASVs that were unique (above 1% relative abundance in that treatment only) to it.

Fig. 4: Summary of the composition, diversity and shared ASVs within sample groupings.
figure 4

Similarity of the composition of microbial communities on different substrata in different environments. Samples are grouped by weighted uniFrac distance using ward linkage (dendrogram) between sample types (colored by environment) and mean community composition at the phylum level is shown. Those phyla that are grouped into ‘Other’ are phyla that are <1% mean relative abundance. Mean relative abundance of several orders that have previously been suggested to be associated with plastics and Simpsons Index of Diversity (showing median and interquartile range) for each group is also shown. The number of ASVs that are shared between different substratum types (that are >1% in relative abundance in these samples) is shown at the bottom.

Full size image

Environmental variables—and not plastic type—have the largest impact on microbial composition

To determine which metadata categories have the largest impact on microbial composition, we constructed random forest models for all 20 categories—including factors such as environment, geographic location, temperature, primer pair, and plastic type to varying degrees of specificity—for each taxonomic level: phylum, class, order, family, genus, species, and ASV (i.e., 140 random forest models were constructed on the rarefied data). We found that models constructed using the microbial composition at the ASV and phylum levels were on average the most (67%) and least (61%) accurate, respectively, and that the light regime (whether samples were incubated at ambient or modified lighting conditions, e.g., shaded or a laboratory 16:8 light:dark cycle) was the metadata category with the highest classification accuracy (maximum 94% at the class level; Fig. 5A). The other metadata categories that were successful more than 80% of the time were: (i) whether the samples came from experiments carried out in the laboratory or the field (maximum accuracy 91% at order or genera level); (ii) whether the samples were incubated/collected from sediment or the water column (maximum accuracy 90% at the ASV level); (iii) the environment that the sample came from (maximum accuracy 86% at order, genera or species level); (iv) the primer pair used for sequencing (maximum accuracy 85% at order or species levels); (v) whether the sample was collected (i.e., unknown environmental residence time) or incubated for a known length of time (maximum accuracy 82% at order, genera or ASV level); and (vi) the DNA extraction method used (maximum accuracy 81% at order level). Depth consistently produced the models with the lowest classification accuracy (maximum accuracy 21% at the ASV level), with temperature also performing poorly (maximum accuracy 40% at the ASV level) and all other categories, such as plastic type, specific incubation time, and geographic location, having intermediate accuracy levels.

Fig. 5: Classification accuracy for random forest models.
figure 5

Classification accuracy (%) for random forest models (classification or regression for discrete or continuous categories, respectively) constructed for (A) all samples grouped within different metadata categories or (B) samples within each environment grouped within plastic type (general) at different taxonomic levels. Random forest models are trained using a subset of 80% of samples (chosen randomly) and classification accuracy is based on testing using the remaining 20% of samples. Figure S4 shows the top most important features at the ASV level across all metadata categories while Supplementary Section 3 shows all taxonomic levels as well as metadata categories.

Full size image

When we examine the taxa with the highest mean feature importance values (proportion that the classification accuracy decreases without that feature), we find that at all taxonomic levels besides ASV, it is a member of the Bacteroidetes with the highest values (maximum 0.159 at the phylum level; Fig. S4 and Supplementary Section 3); even though the Bacteroidetes are substantially less abundant than the Proteobacteria in all but the terrestrial environment (Fig. 4). The next highest feature importance values are generally from the Proteobacteria, and in particular the Alphaproteobacteria (e.g., the ASV with the highest mean value, ASV1197 Erythrobacter, 0.03), although the latitude and longitude usually have the highest individual importance values for both taxonomic groups. There are several metadata categories for which the taxa with high—or low—feature importance values are strongly correlated, for example: latitude and longitude (e.g., ASV1954 Lentimonas and ASV1197 Erythrobacter); study, primer pair, and DNA extraction method (e.g., ASV0715 AEGEAN-169 marine group and ASV0496 SAR116 clade); and source, material type, plastic type and whether the sample comes from an incubation or collection experiment (e.g., ASV0808 Sphingomonadaceae, and ASV4333 Alteromonas; Fig. S4).

The Plastisphere includes potential plastic biodegraders and potentially pathogenic species

To determine the taxa that were potentially specific to plastics, samples were grouped by the environment they originated from and random forest models were trained on general plastic type (as in Fig. 1D). Random forest classification models showed that splitting samples by environment improved the classification accuracy of the plastic type metadata category (to a maximum of 83% at the family, genus or ASV level in the marine environment) and the accuracy remained at below 50% at all taxonomic levels in only the aquatic environment (Fig. 5B and Supplementary Section 3). We were able to identify taxa that were significantly differentially abundant between substratum types in all environments, although in the terrestrial environment we only found significant differences at the genus level (and not at any other taxonomic level; Fig. S5 and Supplementary Section 3). In the marine environment, these included large numbers of Alphaproteobacteria, that were more abundant in planktonic samples, while the Bacteroidetes and Gammaproteobacteria were more abundant in biofilm samples (either plastics or controls). Of particular interest were the taxonomic groups contained within the hydrocarbonoclastic Oceanospirillales (the families Saccharospirillaceae and Halomonadaceae and genera Alcanivorax and Oleiphilus) and Alteromonadales (the families Alteromonadaceae, Marinobacteraceae and Pseudoalteromonadaceae) [118] orders that were more abundant in plastic samples than control biofilms, with the Oceanospirillales generally being more abundant in the aliphatic plastic samples and the Alteromonadales generally being more abundant in the other plastic samples (Fig. S5 and Supplementary Section 3). The Cyanobacteria and the hydrocarbonoclastic Halomonadaceae were always more abundant in the unknown plastic samples (Fig. S5 and Supplementary Section 3), which were collected from the ocean after unknown residence times (Supplementary Table S2).

In the aquatic environment, there were also several potentially hydrocarbonoclastic taxonomic groups that were more abundant in plastic than other samples (the families Thalassospiraceae, Alteromonadaceae, Pseudoalteromonadaceae, Saccharospirillaceae and Xanthomonadaceae and genera Idiomarina and Alcanivorax), however, this difference was only significant for the potentially PAH-degrading Spirosomaceae family [119]. In the freshwater environment, there were also several taxonomic groups that have either been isolated from hydrocarbon-contaminated environments (ASV0388 Arcobacter and ASV0394 Arcobacter cryaerophilus [120]; ASV3841 Unclassified Xanthomonadaceae [118]) or have been suggested to be capable of hydrocarbon (Unclassified Immundisolibacteraceae genera [121]) or biodegradable plastic (ASV2841 Flavobacterium [122]) degradation that were significantly more abundant in the unknown plastic samples. Those that were more abundant in other plastic samples, were classes (Acidobacteriia and Thermoanaerobaculia) and orders (OPB56—Ignavibacteria—Caulobacterales, PB19—Deltaproteobacteria—and Verrucomicrobiales), and it is therefore much more difficult to speculate on the metabolic potential of these taxa. The same was true of the terrestrial environment, where many of the genera that were more abundant in the plastic samples remained unclassified at the order level (e.g., Frankiales and Microtrichales), although several genera known to be able to degrade hydrocarbons were also identified, e.g., Arthrobacter, Acinetobacter [118], Methylocaldum [123], and Nitrosomonas [124] as being significantly more abundant in the plastic (unknown plastics only) than control biofilm samples.

Previous studies have indicated the presence of potentially pathogenic hitchhikers in the Plastisphere, e.g., [30, 64, 125,126,127,128] however, other studies have found that potentially pathogenic species were actually higher in abundance on natural substrates, such as wood, e.g., [25]. Therefore, we investigated taxa that could potentially be pathogenic and that were more abundant on plastics than other samples. We found several taxonomic groups that: (i) were potential pathogens of animals—Tenacibaculum (abundant in marine aliphatic plastic samples) and unclassified Pirellulaceae (abundant in marine other plastic samples); (ii) contained lower taxonomic levels with human pathogens—Clostridiales and ASV0589 Thalassospira (both abundant in marine other plastic samples) and Chlamydiae (abundant in marine unknown plastic samples). However, the ability of a bacterium to be pathogenic is likely dependent on the presence of specific virulence factors which are often in mobile genetic elements [129]. The resolution of amplicon sequencing data is therefore insufficient to determine pathogenesis. There were also taxa, for example the genus Vibrio, of which members are potentially capable of degrading plastics [130] as well as being pathogens of humans and other organisms [129]. Curiously, we found that the Vibrionales were more abundant in plastic samples than control biofilm or planktonic samples in the marine environment, but in the aquatic environment they were more abundant in planktonic than any other samples. This highlights the need for other methods that will enable differentiation between pathogens and non-pathogens, as well as between strains that are capable of plastic degradation and those that are not (as discussed in detail in [131]).

Microbial community succession on different material types

When we examined the taxa that discriminate between early and late incubation times (up to or above 7 days of incubation, respectively), we find some fairly consistent patterns across the marine and aquatic environments: (i) the Bacteroidetes were always significantly more abundant at later time points (where they differed in abundance between early and late time points); and (ii) the Alphaproteobacteria (or orders within the Alphaproteobacteria) are always more abundant at earlier time points (Fig. 6). Other taxonomic groups were more dependent upon the specific comparison. For example, the Gammaproteobacteria, that were more abundant at later time points in marine plastic samples (aliphatic and other plastics), but were generally more abundant at early time points in the aquatic environment, with the hydrocarbonoclastic Oceanospirillales and Alteromonadales in particular being more abundant at early time points. There were also a large number of phyla that were significantly (p < 0.05; Wilcoxon rank sum tests with holm-bonferroni false discovery rate correction) more abundant at early than late time points in the aquatic control biofilm samples, namely the Chloroflexi, Verrucomicrobia, Actinobacteria, Cyanobacteria, and Deltaproteobacteria (Fig. 6). Both other plastic and control biofilm samples in the freshwater environment showed different successional patterns than in either the marine or aquatic environments, although some members of the Bacteroidia were still more abundant at later time points. In the plastic samples, only the Methylophilaceae were more abundant at early than late incubation stages, while in the control biofilm samples there were no taxa that were significantly more abundant at late incubation stages.

Fig. 6: Differential abundance of taxa between early and late incubation times.
figure 6

Heat trees showing differential abundance of taxa between early and late incubation times (up to or above 7 days, respectively) within substratum types (aliphatic, other plastic, and control biofilms) for each of the marine, aquatic and freshwater environments where samples were taken at different time points. Strength of colors indicates differential abundance, with gray indicating no significant difference (p > 0.05; Wilcoxon rank sum tests with holm-bonferroni false discovery rate correction), and strong yellow or red colors indicating that log2 fold change is at least threefold higher in the early or late samples, respectively. Some key taxa are indicated in the empty, larger tree on the left and all are shown in Supplementary Section 4. Tests between different substrata at the same time point as well as between different substrata at all time points are also shown in Supplementary Sections 4 and 5. All terrestrial samples were collected after an unknown environmental residence time and could therefore not be included.

Full size image

The differences between substratum types were also compared at different time points, revealing different patterns between the marine and aquatic environments; there were more differences between sample pairings at late time points in the marine and at early time points in the aquatic environment, with a greater number of differences overall in the aquatic environment (Supplementary Sections 4 and 5). These again included the hydrocarbonoclastic Oceanospirillales and Alteromonadales, that were more abundant in aliphatic plastic samples than control biofilms at early time points, but these differences had largely disappeared by later time points. In the freshwater environment there were not enough samples for comparisons to be made for most sample pairings (only the other plastic and control biofilm samples had both early and late time points), however, differences between sample types were only observed at late incubation times (between other and unknown plastics as well as other plastics and control biofilms).

Discussion

The re-analysis of 35 Plastisphere amplicon sequencing studies that we present here has allowed for a comprehensive characterization of the 16S rRNA gene community that colonizes plastics in marine, freshwater, other aquatic and terrestrial environments. We trained random forest models to classify samples within metadata categories, revealing that, overall, there were a number of variables—both relating to study design and environmental factors—that were better able to differentiate between sample groupings than plastic type (Fig. 5A). This agrees with the small re-analysis of studies presented by Oberbeckmann and Labrenz (2020; n = 5) [45], who reported that geographic location and salinity, but not substrate type, significantly discriminated between different samples. They did not, however, investigate different aspects of study design, such as environment, light regime or DNA extraction methods and only included studies that used the same primer pair, all of which are categories that we found to have a higher classification accuracy than location (latitude and longitude) or salinity (Fig. 5A). We tried to control for primer pair and DNA extraction technique by agglomerating sequences based on phylogenetic placement, which is known to reduce the variation between different 16S rRNA gene regions [50], although there are some taxa that will not be amplified by particular primer pairs, including the hydrocarbonoclastic groups targeted by many Plastisphere studies [132].

When we split our analyses to different environments (marine, aquatic, freshwater, and terrestrial), we find that plastic type could be used to train random forest models that were up to 83% accurate (Fig. 5B). We were, therefore, able to show that the differences identified by individual studies between microbial biofilms on control substrates and the Plastisphere were the same across the Plastisphere studies that we re-analysed here, regardless of the methodological differences between studies. For example, we find that the hydrocarbonoclastic Oceanospirillales and Alteromonadales, both previously suggested to be plastic biodegraders, e.g., [35, 73], were both important for differentiating between different substrate types in the random forest models (i.e., between aliphatic plastics, other plastics, control biofilms, and planktonic communities). They were also identified as significantly differentially abundant between these substrates. Furthermore, we find that the Oceanospirillales are higher in abundance in aliphatic plastics while the Alteromonadales are higher in abundance in the other plastic types. This is an observation that, to our knowledge, has not been previously acknowledged. Members of both of these orders, e.g., Alcanivorax sp. (Oceanospirillales) [35] and Alteromonas sp. (Alteromonadales) [133], have been reported to be able to degrade PE (aliphatic plastic) or PET (other plastic), respectively. However, these studies do not incorporate measurements of the plastic carbon assimilation into biomass (or of respiration) and therefore fall short of being definitive measurements of plastic biodegradation [134].

Almost 1,000 toxic chemicals are known to be associated with plastics [135]. These include manufacturing additives as well as weathering sub-products, most of which are not chemically bound and leach from plastics on environmental exposure [136]. They are more labile and are therefore more likely to be biodegraded by the Plastisphere than plastics themselves. However, there are only two studies included here that consider the effect that these have on community composition [31, 33]. For many plastics, these chemicals are not only additives but are also likely intermediates of both biotic and abiotic degradation. Cyclic and noncyclic dicarboxylic acids, for example, are the most commonly identified intermediates for the abiotic degradation of PE, PP, PS, and PET [137] and are also intermediates of biotic PET [138] and polyester [48] as well as plastic additive (plasticizer) [37] degradation. It has been suggested that the microbial community that is likely able to use these more labile chemicals is only present at earlier time points [31]. Almost all surfaces that enter the environment go through distinct and well-characterized stages of microbial community succession, and it is well known that microbial communities on different surfaces converge at later time points [139,140,141]. Here we have defined early incubation stages as before and later as after 1 week of incubation. However, due the relatively small number of studies that characterize the Plastisphere across a wide range of different incubation times, we do not actually know what early or late incubation stages are within the Plastisphere. Long residence times in the order of weeks to months have generally been favored by Plastisphere studies—the mean incubation time for all samples included here was 100 days—rather than the days that it typically takes for microbial communities to diverge from organisms that are efficient at degrading the surfaces they colonize to cheaters, cross-feeders, and grazers [139, 142]. This makes sense when we know that plastics may take longer than hundreds of years to be completely mineralized [143], but it also explains why several studies have found that potential degraders are only present at very low relative abundances [26] or are abundant only at early time points [31].

There are also several other areas that we have identified as needing further research before definitive conclusions may be made. For instance, there are currently only four studies performed in the terrestrial environment (and only two that could be included), all current Plastisphere studies were performed in the Northern Hemisphere and the majority were performed in temperate environments (Fig. 1). Whilst this hinders our ability to draw large-scale conclusions, it offers an opportunity for researchers to ensure that future data collected in these areas are comparable and address these knowledge gaps. We were unable to confirm or reject the suggestions that plastics carry higher abundances of pathogenic hitchhikers, e.g., [23, 30, 64, 126,127,128], and/or antimicrobial resistance genes, e.g., [39], than control biofilms. We do note, however, that due to their recalcitrance and buoyancy plastics present a higher chance of carrying pathogens or antimicrobial resistance genes across greater distances than most natural surfaces. We, along with many of the studies included here, have identified the presence of potential plastic biodegraders and potentially pathogenic organisms on plastics. However, the amplicon sequencing used does not usually give the resolution required to differentiate between closely related strains of the same species [144], for example between pathogenic and non-pathogenic Vibrio spp. [129, 145] or hydrocarbon degrading and nondegrading Pseudomonas putida (previously P. oleovorans) spp. [146]. The presence of genes encoding the production of enzymes used in hydrocarbon biodegradation [118], or specific virulence factors required for pathogenicity are often conferred by mobile genetic elements [129] and their presence is therefore not necessarily based in phylogeny. The use of metagenomics rather than amplicon sequencing in future studies would aid in determining whether these potential plastic biodegraders and potential pathogens are indeed biodegraders or pathogenic [131]. As discussed further in [131], it may even generate candidate enzymes that could be tested for their ability to degrade plastics, such as in Danso et al. [147] for PET hydrolases.

This study presents, for the first time, a comprehensive re-analysis of all Plastisphere studies that utilize amplicon sequencing of the 16S rRNA gene. We have revealed through machine learning methods that environmental factors, such as environment and light availability as well as aspects of study design, such as primer pair and incubation time play a large role in shaping Plastisphere community composition. Notably, we have identified members of the microbial community that are consistently more abundant in biofilms formed on plastics than control biofilms across multiple studies and environments. This highlights the urgent need to determine whether these microbes are capable of plastic biodegradation or are pathogens of humans or other organisms. We have also identified a number of other key areas in which we are lacking even basic knowledge and where future research should be directed. It is clear that plastic pollution is a key indicator of the Anthropocene and we must focus future research on gaps in our knowledge that we highlight here if we wish to mitigate its effects.