A Comparative Taxonomic Profile of Microbial Polyethylene and Hydrocarbon-Degrading Communities in Diverse Environments.

BACKGROUND: Polyethylene (PE) is one of the most abundant plastic wastes which accumulates in marine and terrestrial environments. As microbial degradation has been a promising approach for the bioremediation of polluted environments, identification of the microbial community profile where these pollutants accumulate, has recently been in focus.
OBJECTIVE: We have investigated the taxonomic and functional characteristics of polyethylene- degrading microorganisms in a plastic waste recycling site in Tehran, Iran.
MATERIALS AND METHODS: We have analyzed and compared a 16S rRNA dataset from this study with 15 datasets from 4 diverse plastic and oil polluted habitats to identify and evaluate bacterial communities involved in bioremediation.
RESULTS: Our findings reveal that Proteobacteria, Actinobacteria, Acidobacteria and Cloroflexi were the dominant phyla and Actinobacteria, Alphaproteobacteria, Gammaproteobacteria and Acidimicrobia were dominant classes in these samples. The most dominant Kegg Orthology associated with PE bioremediation in these samples are related to peroxidases, alcohol dehydrogenases, monooxygenases and dioxygenases.
CONCLUSIONS: Long-term presence of contaminants in soil could lead to changes in bacterial phyla abundance, resulting in metabolic adaptations to optimize biological activity and waste management in a diverse group of bacteria. Copyright:

1. Background

One of the challenges of our planet in recent years has been the accumulation of plastic debris in terrestrial and marine ecosystems through anthropogenic activities leading to extensive contamination of the environment ( 1 ). Plastic is a broad name for different types of synthetic polymer with high molecular weight and is generally derived from various petrochemical and hydrocarbon sources ( 2 , 3 ).

Since polyethylene (PE) consists of long chains of hydrocarbon with very balanced charge and various degrees of branching and also lacks chemical groups ( 4 ), it has become a stable polymer with widespread uses in different industries and agriculture ( 5 ). The most consumed synthetic polymer is PE with a current global production of 140 million tons per annum in 2001 ( 6 ). Global demand for polyethylene resins was 99.6 million metric tons in 2018 (www.freedoniagroup.com).

Since 2015, approximately 8.3 billion metric tons of plastic have been produced out of which 6.3 billion metric tons have become waste. Of this waste 12% was incinerated, 79% was stored in landfills and only 9% has been recycled ( 7 ).

Microorganisms have important and key roles in their various ecosystems ( 8 ) and are involved in the degradation of natural and man-made polymers ( 9 ). Standard culturing techniques can be used to isolate and identify only less than 1% of microorganisms in many environments. Therefore, culture-independent techniques are now a better tool in the detection of millions of unknown species in different environments with applications in biotechnology ( 10 ) using molecular ecological methods to overcome this limitation ( 11 ).

Marker gene amplification metagenomics using PCR amplification and sequencing marker genes such as the 16S rRNA genes, provide extensive information about community and taxonomic distribution profile of microorganisms and their probable functions ( 12 - 14 ).

2. Objectives

In this study, taxonomic characteristics of polyethylene- degrading microorganisms in a plastic waste recycling site in Iran was investigated. A 16S rRNA dataset from this study as well as 3 contaminated petroleum sites and 12 samples of PE contaminated soils from SRA databases for comparison were used.

Materials and Methods

3.1. Sample Collection and Enrichment

Soil samples were collected by random sampling method ( 15 ) from plastic waste recycling site at a Tehran city suburb in Iran (35ᵒ30’22.1”N 51ᵒ22’51.3”E) ( 16 ) with a semi-arid climate. The samples were transferred to the laboratory on ice and processed within 24 hours.

The enrichment procedure was carried out to assist in the survival of microorganisms capable of consuming and degrading PE as the sole carbon source. Therefore, 10 g of each soil samples was placed in a flask containing 20 mL Bushnell Hass (BH) broth medium (per liter of distilled water): 0.2 g MgSO4, 0.02 g CaCl2, 1 g KH2PO4, 1 g K2HPO4, 1 g NH4NO3, 0.05 g FeCl3 with a final pH of 7. Four flasks for each sample contained, as the sole carbon source, 1% LDPE powder (Low Density Polyethylene), 1% LDPE powder and 0.1% (yeast extract + glucose), 1% hexadecane (as model substrate) ( 17 , 18 ) or 1% hexadecane and 0.1% (yeast extract + glucose). The flasks were incubated in a shaker at 37ºC, with an rpm of 120 for 30 days and subsequently mixed for metagenomic DNA extraction.

Fifteen 16S rRNA datasets publicly available in SRA (Sequence Read Archive) from 3 different environments were used as control and for comparison to the sample from this study. Six samples (ISP; Italy Soil PE) were collected at dump site in Fiorenzuola d’Arda in Northern Italy with temperate climate and average temperature of 12.8ºC and also mean annual rainfall of 774 mm ( 19 ) . Six samples (CSP; China Soil PE) were collected from cotton field in Shihenzi city with a desert climate and an average temperature of 7.8ºC and also mean annual rainfall of 225 mm ( 20 ). Three samples (PSO; Poland Soil Oil well) were collected at an oil well site in Weglowka near Krosno Poland with cold and temperate climate with the average temperature of 6.9 ᵒC and mean annual rainfall of 816 mm ( 21 ). Our sample (TSP; Tehran Soil PE) was collected in a waste recycling workshop at Kahrizak suburb in Tehran with cold-desert climate and average temperature of 16.9 ºC and 156 mm average annual rainfall (https://en.climate-data.org/).

3.2. Total DNA Extraction and Illumina Sequencing

In order to extract total metagenomic DNA from soil samples, the following extraction protocol was used: 4 flasks of each sample were mixed following which 1 g sludge was added to 200 µL extraction buffer (NaCl 1.5 M, Tris-HCl 100 mM, EDTA 100 mM, pH=8), 30 µL NaH2PO4 (0.12 M, pH=7.4) ( 22 ), 0.1 g sterile glass powder ( 23 ), 30 µL SDS 25% and 20 µL Tween 20. The sample was vortexed for 5 min and incubated at 65ºC for 15 min (with invert mixing after every 1-2 min). The tubes were centrifuged at 15000 x g for 10 min at 4ºC to collect the supernatants which were then mixed with 10 mg.mL -1 lysozyme, 20 mg.mL-1 proteinase K followed by incubation at 37ºC for 1 h and inverting the tubes every 10 min. In the next step, equal volume of chloroform: isoamylalcohol (24:1 v/v) was added to every tube and centrifuged at 12000 x g for 15 min at 4ºC to collect the upper phase. DNA was precipitated by adding ice-cold absolute ethanol and centrifugation at 2000 g for 20 min at 4ºC. Supernatant was discarded and the pellet was washed with 70% ethanol and centrifuged at 2000 g for 20 min at 4ºC. Finally, the pellet was air-dried and dissolved in 50 µl TE buffer. Concentration and quality of the extracted DNA were determined using NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, USA) and 1% agarose gel electrophoresis, respectively.

The V3-V4 variable regions of 16S rRNA genes were amplified and sequenced by Illumina Miseq Paired-end Technology (Macrogen, Korea). The raw sequence reads were deposited in to the NCBI sequencing read archive (SRA) under Accession No PRJNA597282.

3.3. Data Analysis

Raw sequences were analyzed using QIIME2-2019.7 ( 24 ). The sequence reads were demultiplexed and denoised with DADA2 ( 25 ) based on quality scores and filter

3.4. Functional Metagenomic Prediction Using 16S rRNA Datasets

Prediction of functional profile was carried out using PICRUSt ( 27 ). After OTU data generation for all 16S rRNA data based on the GREENGENES reference database, biome file format was prepared as input for PICRUSt2 software.

Full pipeline for EC number, KEGG Orthology (KO) and MetaCyc pathway abundance prediction was carried out using picrust2_pipeline.py script. The Nearest Sequenced Taxon Index (NSTI) was assumed >2 by default. Phylogenetic placement of reads was done by HMMER, EPA-NG ( 28 ) and Gappa ( 14 ) algorithms, KO prediction by Castor ( 29 ) algorithm and relative pathway abundance by MinPath ( 30 ) implemented in HUMAnN2 ( 31 ).

4. Results

4.1 Bacterial Community Structure

For a comprehensive study of the microbial community that is capable of degradation and PE consumption, enriched samples collected from Tehran were compared with 12 samples of PE contamination at a cotton farm in China, landfill in Italy and 3 samples of oil well soil in Poland that were previously published and deposited in NCBI, SRA (Table 1).

Table1

Summary of 16S rDNA datasets used for analysis

SRA Accession Number	ID*	Treatment	Location	Duration of Exposure	Total Reads	High quality reads	Reference
SRX7433929	TSP	PE Enrichment	Iran	20 years<	97543	13807	This study
SRX4084527	PSO1	Oil Well	Poland	30 years<	49336	10883	( 20 )
SRX4084525	PSO2	Oil Well	Poland		43910	8346
SRX4084523	PSO3	Oil Well	Poland		52014	11575
SRX5439652	CSP1	Soil Around PE	China	30 years<	67144	17678	( 19 )
SRX5439651	CSP2	Soil Around PE	China		49933	13291
SRX5439646	CSP3	Macro PE Enrichment	China		40122	10171
SRX5439645	CSP4	Macro PE Enrichment	China		32687	6349
SRX5439643	CSP5	Micro PE Enrichment	China		65432	18082
SRX5439642	CSP6	Micro PE Enrichment	China		34382	4462
SRX6798559	ISP1	PE Plastic	Italy	35 years<	23678	5791	( 18 )
SRX6798558	ISP2	PE plastic	Italy		28467	7391
SRX6798553	ISP3	Control Soil	Italy		17803	5846
SRX6798550	ISP4	Control Soil	Italy		13422	4423
SRX6798549	ISP5	Dump Soil	Italy		14674	3842
SRX6798548	ISP6	Dump Soil	Italy		14803	4099

The IDs used in the table are not taken from the relevant articles reference and are as follows: TSP (Tehran Soil PE sample), PSO (Poland Soil Oil well samples), CSP (China Soil PE samples), ISP (Italy Soil PE

Data analysis were performed using Qiime2 software, following which 4202 OTU were obtained consisting of 36 phyla, 90 classes, 140 orders, 184 families, 200 genera and 50 species. Alpha diversity in rarefaction curves based on observed OTU showed acceptable and nearly complete sampling for all of the samples (Fig. 1).

Figure 1

Rarefaction curve based on observed OTU for 16 samples show acceptable sampling with median frequency 7800 as maximum sampling depth

Bacterial community profile of samples is shown in the phylum and class level in bar plot (Fig. 2A and B). Dominant phyla based on the mean relative abundance include Proteobacteria 27.5% (9.5-49.11%), Actinobacteria 27.24% (7.28-67.60%), Chloroflexi 9.039% (0.829-17.75%) and Acidobacteria 6% (0.511-23.091). The main classes in all samples are Actinobacteria 21.5% (2.7-66%), Alphaproteobacteria 19.45% (8.9-40.4%), Gammaproteobacteria 7.4% (0.06-27.66%) and Acidobacteria 3% (0.2-8.1%) mean percentage abundance.

Figure 2

Taxonomic distribution of microbial communities A At the phylum level. B At the class level (the classification with less than 1% categorized in "others")

In addition, Micrococcaceae, Nocardiodaceae, Streptomycetacease and Mycobacteriaceae in the Actinobacteria class and Sphingomonadaceae, Rhodospirillaceae and Rhodobacteraceae in the Alphaproteobacteria class were the prevailing families. The Gammaproteibacteria class was dominant with Alcanivoracaceae, Pseudomonadaceae and Xanthomonadaceae families in TSP and PSO samples.

The bacterial community is mostly represented by members of the genera Mycobacterium sp., Streptomyces sp., Rhodococcus sp., Nocardia sp., Microbispora sp. and Cellulomonas sp. in the Actinobacteria class and Kaistobacter sp., Sphingomonas sp., Skermanella sp., Rubellimicrobiom sp., Hyphomicrobium sp., Devosia sp. and Agrobacterium sp. in the Alphaproteobacteria class, Pseudomonas sp., Alcanivorax sp., Pseudidiomarina sp., Crenothrix sp. and Thermomonas sp. in the Gammaproteobacteria class. Most of the mentioned genera are active in petroleum and n-alkane degradation.

Differences in abundance and diversity of taxonomies in various samples are dependent on many different factors. Bacterial profiles in various environment follow different patterns which may be related to environment temperature, latitude, nutrient availability, soil structure and the ecosystems ( 21 , 32 , 33 ). Beta diversity results are specified on heatmap for frequencies above 50 abundance (Fig. 3A) and three dimentional Principal Coordinates Analysis (PCoA) plot based on weighted and unweighted unifrac of community profiles based on different treatments in diverse samples (Fig. 3B).

Figure 3

Beta diversity. A The heatmap of the relative abundance of bacteria with frequency above 50 in diverse samples (Sample name are in rows and taxonomies at class level are on columns). B PCoA plot of different treatment in diverse samples, B1 weighted-unifrac and, B2 unweighted-unifrac.

4.2. Prediction of Functional Profile

To evaluate the functional composition and potential of the microbial community, sequence data from 16S rRNA genes of all samples was used in PICRUSt software to predict metabolic structures. Predicted enzymes and proteins were classified by KEGG Orthology (KO) ( 34 ). 7610 KO and 428 metabolic pathways reconstructed by HUMAnN2, were implemented in PICRUSt. 428 pathways were identified in MetaCyc analysis across all samples, among which 48 pathways are involved in degradation of aromatic, aliphatic and xenobiotic compounds.

These 428 pathways were categorized in 14 pathway classes: 1) cofactor, vitamin, electron carrier, prosthetic group biosynthesis, 2) secondary metabolite metabolism, 3) generation of precursor metabolite and energy, 4) amino acid metabolism (biosynthesis and degradation), 5) C1 compound utilization and assimilation, 6) carbohydrate metabolism, 7) nucleoside and nucleotide metabolism (biosynthesis and degradation), 8) fatty acid and lipid metabolism (biosynthesis and degradation), 9) amine and polyamine metabolism, 10) cell structure biosynthesis, 11) inorganic nutrient metabolism and 12) aromatic, aliphatic and xenobiotic degradation, 13) nucleic acid processing, 14) other biosynthesis. Most pathways are involved in vital activities of bacteria and a low percentage of predicted pathways are involved in bioremediation if necessary (Fig. 4A).

Figure 4

Predicted functional profile. A Pie chart showing mean percentage of metabolic pathway classes of microbial communities. B Similarity between predicted function of different treatment sample based on Bray-Curtis dissimilarity distance.

Predicted metabolic profiles of samples show most pathways and gene families are involved in vital cell activities ( 33 ) which is not surprising since microorganisms try to adapt to the fluctuations and stresses in the environment. Therefore, they must tune their metabolic pathways to allow for the consumption and degradation of these pollutants to survive in these environments. Bray-Curtis dissimilarity distance based on predicted KOs and MetaCyc pathway between different samples (treatments) are shown in PCoA plot (Fig. 4B).

In this study, relative abundance percentage of mentioned predicted KEGG Orthology in various samples were calculated (Fig. 5).

Figure 5

Relative aboundance percentage of notable KO in alkane and aliphatic degradation in various samples.

5. Discussion

Microbial degradation of xenobiotic waste from the environment is a promising approach in bioremediation. However, identification and characterization of the microbiota profile in the polluted environment is a priority and a necessity ( 33 ). As the bioinformatic analysis of the microbial community profiles from different soil samples and PE biofilm specimen elucidate, Proteobacteria, Actinobacteria, Chloroflexi and Acidobacteria are the dominant phyla in the compared samples; these most abundant bacteria have also been reported in studies on long term diesel contaminated soil ( 35 , 36 ).

Proteobacteria was the most dominant phylum in control (ISP 3-4) and enriched soil samples (TSP, CSP 1-2, PSOs) but in PE biofilm specimens, Actinobacteria were predominant. According to research of Janssen and coworkers (2006), as expected, the most abundant phylum in soil are Proteobacteria ( 37 ). Militon and colleagues (2010) showed that before and after a 12 month period soil treatment with hydrocarbon, the most abundant phylum in soil were still Proteobacteria, but after a 24 month period, Actinobacteria became dominant ( 38 ). These results confirm that with changes in available nutrients over time, soil bacterial community structure changes and Proteobacteria is typically the dominant phylum in soil. Moreover, it seems that Actinobacteria play a key role in biofilm formation and bioremediation of xenobiotic compounds.

In addition to the mentioned dominant taxa, bacterial communities are diverse in all samples and according to Bell and coworkers (2005), in order to complete metabolic cycles, diverse communities prepare and perform better in ecosystems ( 39 ). A study by Delgado-Baquerizo and colleagues has, for the first time, identified the most abundant bacteria with a collection of soil samples from around the world and suggested these dominant taxa account only for 2% of soil diversity and this small portion of the total diversity forms a major fraction of the soil sample sequences worldwide ( 40 ). Therefore, some taxa are present and abundant in soil samples ; however, a large fraction of diverse taxa present in soil are discounted for even though they may play pivotal and important roles in environment ( 41 ).

PE is an n-alkane polymer with petroleum origin and consists of high molecular weight hydrocarbon with C-C and C-H bonds ( 4 ). Many microorganisms can degrade alkanes as energy and carbon source. The most important enzymes involved in the degradation of environment pollutants, encoded by major gene families, are methane monooxygenases, soluble cytochrome P450s, non-heme iron monooxygenases (alkane hydroxylase, alkB), alcohol and aldehyde dehydrogenases and dioxygenases ( 42 , 43 ). Ligninolytic enzymes such manganese peroxidase and laccase can cause oxidation and molecular weight reduction of PE ( 4 , 44 - 48 ). All of mentioned enzymes belong to the family of oxidoreductases; these enzyme families, with different electron acceptors, have notable potential for degradation of the xenobiotic compounds ( 49 ). The most abundant KO relative with PE degradation in all samples belong to peroxidases, alcohol dehydrogenases, dioxygenases and monooxygenases.

6. Conclusion

The accumulation of plastic pollutants in terrestrial and marine ecosystems is a major environmental problem. Many microorganisms with ability to degrade and eliminate these contaminants have been identified by available culture-dependent methods. Metagenomic approaches provide extensive information on the microorganism populations, their relationship with the environment and can consequently, identify their roles in the process of polyethylene degradation. In this study, we re-identified some of active bacterial genera in n-alkanes bioremediation such as Mycobacterium sp., Rhodococcus sp., Streptomyces sp., Nocardia sp., Sphingomonas sp., Agrobacterium sp., Alcanivorax sp. and Pseudomonas sp. The most dominant KO in these samples are related to peroxidase, alcohol dehydrogenase and dioxygenase pathways that exist in the mentioned phyla. A deeper understanding of the complex interactions between different bacterial communities and environmental milieu can expedite and improve extant bioremediation processes.