Lignocellulose-mediated selection of potential halophilic PET-degrading enzymes from mangrove soil

Domain: chemistry
Published: 2026-04-07
Editorial outcome: minor_revision

Abstract
Mangroves are ecosystems located at land–sea transition zones, where they are continuously exposed to plant biomass and plastic pollution. Their soils harbor extensive microbial diversity with potential for discovering polymer-degrading enzymes. Here, we perform a microcosm experiment to examine how mangrove soil microbial communities respond to inputs of lignocellulose or polyethylene terephthalate (PET) in the presence and absence of seawater, and to explore the selection of putative PET-active enzymes (PETases) using gene- and genome-resolved metagenomics. Incubation conditions lead to a gradual increase in salinity, resulting in the enrichment of halophilic taxa, including spore-forming bacteria and archaeal species, particularly in seawater-depleted treatments. Lignocellulose input is the primary driver of soil microbial community restructuring, followed by seawater presence. In dry, lignocellulose-amended microcosms (L treatment), microbial diversity is significantly reduced, while lignocellulolytic taxa within the phyla Bacillota and Actinomycetota are enriched. Twelve potential PETases are identified in the L treatment, sharing >70% sequence similarity with known PETases, and three are predicted to be thermostable. Two putative PETases from Microbulbifer species display distinct sequence and structural features, thereby expanding the currently limited PETase sequence landscape. This study demonstrates that perturbing environmental microbiomes with plant-derived polymers represents a promising strategy for capturing novel PETases.

Introduction
Mangroves are critical ecosystems that provide multiple global services, including carbon sequestration, climate change mitigation, and biodiversity support (1,2). Located at the terrestrial-marine interface, these ecosystems are dynamic and shaped by tidal currents, hosting diverse microbial species that participate in complex interactions, metabolic processes, and ecological functions (3,4,5). The microbial communities inhabiting mangrove soils are predominantly composed of halotolerant species dispersed through seawater floods (6,7,8). Mangrove soils are prone to multiple environmental changes and disturbances. For example, (i) constant flooding introduces new microbial species and promotes nutrient runoff (9,10); (ii) pollutant inputs (e.g., microplastics and agricultural residues) can modify physicochemical conditions and alter microbial community structures (11,12); and (iii) drought events increase soil salinity, which reduces plant productivity and elevates mortality rates (13,14).

Globally, plastic pollution represents a major environmental challenge, with polyethylene terephthalate (PET) being a significant contributor to marine waste. Mangrove ecosystems, as natural sinks, are particularly vulnerable to the accumulation of fossil-derived plastics (15,16). Exposure of mangrove soil microbiomes to tidal flooding, salinity shifts, plant-derived polymers (17,18), and microplastics (19,20) makes these ecosystems promising reservoirs of novel lignocellulose- and plastic-transforming enzymes with polyextremophilic traits potentially suitable for industrial applications (21,22). Halophilic enzymes, such as lipases and (hemi)cellulases, are valuable for numerous biotechnological applications, including their use as active compounds in detergents and as components of enzymatic cocktails for the saccharification of agricultural residues in biorefineries (7,23). However, halophilic PET-depolymerizing enzymes from mangrove-associated environments have not yet been described. Recently, we reported five putative PET-depolymerizing enzymes from Pseudoxanthomonas winnipegensis, an abundant member of PET-transforming bacterial consortia selected from a mangrove soil microcosm (24).

PET-depolymerizing enzymes are α/β-fold hydrolases (known as PETases), primarily found in taxa within the phyla Pseudomonadota, Actinomycetota, and Bacillota (25,26). Based on amino acid sequences and 3D structural conformations, these enzymes have been classified into three main groups (I, IIa, and IIb) (27). More recently, a fourth group of marine-derived PETases (group III) has been proposed, encompassing salt-tolerant enzymes from the Halopseudomonas lineage (28). Additionally, three halophilic PETases with high homology to IsPETase (from Piscinibacter sakaiensis, formerly Ideonella sakaiensis) (29), have been identified in protein catalogs from deep-sea sediment metagenomes (30,31). From an industrial perspective, PETases that remain active under extreme conditions, such as high temperatures (~70 °C), low pH, highly crystalline PET, and elevated salt concentrations, are particularly desirable for optimizing biocatalytic processes (26,31,32,33).

Ecological perturbation experiments enable the investigation of dynamic responses and resilience of microbial communities under acute or chronic physicochemical changes (34). However, these experiments with soil microbiomes (e.g., microcosms) can be adapted to select for target or rare microbial populations (35) with PET-depolymerizing capabilities (24,36). For example, a recent two-year microcosm experiment using deep-sea sediments and PET microparticles successfully selected PET-degrading prokaryotic species (37). Although PET microparticles or films are rational inducers of such selection, other easy-to-break PET analog substrates (e.g., plant biomass or polyester-rich substrates) can also serve as effective disturbance factors, increasing the likelihood of enriching PET-depolymerizing taxa and enzymes (38). This is because natural substrates of PETases, such as plant polymers (e.g., xylan, suberin, and cutin), contain the same ester bonds present in PET (39,40). Therefore, selecting for lignocellulose-degrading microbial communities in soil microcosms represents a promising strategy for enriching prospective PET-active enzymes. This idea is further supported by the frequent occurrence of polyester-degrading enzymes in plant polymer-rich environments such as compost, soils, and herbivore gastrointestinal tracts (40,41).

In this work, we design a microcosm experiment to (i) evaluate the responses of mangrove soil microbial communities (i.e., shifts in diversity and structure) to inputs of lignocellulosic substrate or amorphous PET microparticles in a saline environment with or without seawater, and (ii) assess the selection of polymer-degrading prokaryotic taxa and enzymes. By reconstructing metagenome-assembled genomes (MAGs) and generating comprehensive gene catalogs, we screen the enzymatic potential of these restructured soil microbial communities for the transformation of plant biomass and PET. In this study, we demonstrate that ex situ perturbation of mangrove-derived soils with lignocellulose substantially reshapes microbial communities, enriching halophilic species with the capacity to transform plant polymers. The inclusion of lignocellulose in mangrove soil microcosms is particularly associated with the selection of prokaryotic lineages that harbor an expanded genetic potential for PET depolymerization, including putative halophilic PETases.

Results
Reshaping of mangrove soil microbial communities in lignocellulose-amended microcosms
A soil microcosm experiment was conducted with lignocellulose (rice husk) and PET particles amendments, along with two environmental factors: (i) the presence of external seawater-derived microbial populations and nutrients, and (ii) a gradual increase in soil desiccation and salinity over the incubation period (30 and 90 days) at 30 °C (Fig. 1). Based on Shannon index values derived from 16S rRNA gene and ITS2 amplicon sequencing, microcosms containing lignocellulose without seawater (L treatment) exhibited a trend toward lower alpha diversity compared to other treatments without seawater after 30 days of incubation (ANOVA; p > 0.08, Fig. 2A). A significant pairwise difference in bacterial Shannon diversity was observed between the L treatment and controls with seawater (CW), with the latter showing higher diversity at 30 days (p < 0.05). After 90 days, bacterial alpha diversity decreased significantly in microcosms without seawater compared to those with seawater inclusion (p < 0.001, Fig. 2A). Fungal alpha diversity in the L treatment displayed a distinctive pattern, with an initial decline at 30 days followed by an increase at 90 days, suggesting recovery and re-assembly of the fungal community after perturbation.

Beta diversity analysis based on Bray-Curtis distances revealed that both bacterial and fungal communities were significantly reshaped by lignocellulose input (PERMANOVA; p < 0.001, R2 = 0.187 and 0.199), with seawater inclusion also contributing significantly to community variation (p < 0.001, R2 = 0.122 and 0.102, Fig. 2B). These compositional changes were reflected in the relative abundance of dominant phyla after 90 days, including Bacillota (in L and LW), Bacteroidota (in LW), and Planctomycetota (in L) (ANOVA; p < 0.001, R2 = 0.057, Fig. 2C). Furthermore, combined bacterial and fungal co-occurrence network analysis revealed that lignocellulose input substantially reduced the complexity of the soil microbial community after 90 days, as evidenced by a decrease in the number of nodes and edges compared with negative controls (C and CW) and the original mangrove soil microbial community (S) (Suppl. Data 1, Fig. 2D).

Enrichment of spore-forming bacterial taxa in lignocellulose-amended microcosms
The Boruta algorithm identified 102 bacterial ASVs associated with either a specific treatment group or the original mangrove soil bacterial community (S) (Fig. 3A). Among these, taxa within the families Bacillaceae, Sporolactobacillaceae, and Paenibacillaceae showed higher relative abundances in lignocellulose-amended microcosms (L and LW) after 90 days of incubation. Most of these ASVs were enriched in the lignocellulose treatments and exhibited lower abundances in PET-amended microcosms (P and PW), negative controls (C and CW), and in the original soil bacterial community (S) (Fig. 3A). Differential abundance analysis using DESeq2 revealed a significant enrichment (p < 0.0001; log2-fold change > 5) of 86 ASVs in lignocellulose-amended microcosms compared with negative controls (Fig. 3B). Notably, some ASVs assigned to Bacillaceae (3 ASVs), Micromonosporaceae (1 ASV), Paenibacillaceae (1 ASV), and Pirellulaceae (1 ASV) displayed particularly high relative abundances, with fold-changes greater than 9 in lignocellulose-amended treatments over the 30–90 day incubation period.

Differences in abiotic factors modulate microbial species selection
Metagenome sequencing was performed using three biological replicates per treatment (L, LW, P, PW, C, and CW) after 90 days of incubation. This sample selection was supported by detected drastic changes in microbial community structure and diversity after 90 days of incubation (Fig. 2). On average, 65 million reads per sample (ranging from 50 to 150 bp) were obtained, yielding approximately 23.4 Gbp of data per treatment (Suppl. Data 2). Taxonomic classification of metagenomic reads using Kraken2 and Bracken assigned ~27% of high-quality reads to known taxa (Fig. 4AB; Suppl. Data 2). Among the classified reads, 91% corresponded to prokaryotic sequences, with an average of 710 prokaryotic families detected per sample. Bar plot analysis revealed that Streptomycetaceae and Micromonosporaceae were the most abundant families (~20% of assigned reads) in the L treatment, whereas the halophilic archaeal family Haloarculaceae dominated (~15% of assigned reads) in the C treatment. The relative abundance of Haloarculaceae was also higher in the P treatment compared with PW and CW treatments (Fig. 4B), indicating that the absence of seawater and soil desiccation favors halophilic archaea, while lignocellulose addition promotes actinobacterial growth. Shannon index analysis of metagenomic read-based taxonomic profiles showed significantly lower diversity in the L treatment and higher diversity in the LW treatment (p < 0.05, Fig. 4C). Clustering analysis further revealed distinct microbial community structures among the L, LW, and C treatments (Fig. 4D).

Microbial diversity captured by metagenome-assembled genomes
High-quality metagenomic reads from each sample were pooled by treatment and subjected to co-assembly, yielding an average of 2.6 million contigs per treatment (N50 ≈ 952 bp). Approximately 28% of high-quality reads mapped to contigs (Suppl. Data 2). The binning process grouped 141,000 contigs into 196 bins (Fig. 4A). After refinement and dereplication, 66 medium-quality MAGs (≥50% completeness, ≤10% contamination) and 8 high-quality MAGs (≥90% completeness, ≤5% contamination) were retained. In total, ~5% of metagenomic reads mapped back to these 74 MAGs, indicating limited diversity capture. The average MAG size was 3.33 ± 1.34 Mbp, with ~535 contigs per genome (Suppl. Data 3). Taxonomic classification with GTDB-Tk assigned 4 MAGs to Archaea (P1_Haladaptatus, P10_Halomarina, CW2_Nitrososphaeraceae, and CW3_Nitrososphaeraceae), while the remainder were assigned to Bacteria. Interestingly, 5 MAGs enriched in the CW and LW treatments were affiliated with the taxa UBA6522 (an unclassified Gammaproteobacteria), and 4 MAGs enriched in the P and C treatments were affiliated with UBA5704 (an unclassified Acidimicrobiia). Only 18 MAGs were assigned to known and named genera or species in GTDB (Suppl. Data 3).

Clustering of genomes provides insight into their provenance and selection in microcosms
Based on read mapping and relative abundance data, the 74 MAGs were clustered into six significantly different groups or functional guilds (Kruskal–Wallis with Bonferroni correction; p < 0.05) (Fig. 5AB). Group A included MAGs enriched in the lignocellulose-amended treatments (L and LW). Three MAGs (L1_Tuberibacillus calidus, L3_Bacillus smithii, and L2_Cohnella laeviribosi) were particularly abundant in both treatments. Group B comprised MAGs enriched in CW and PW compared with L and LW, including four affiliated with the taxon UBA6522. Group C contained MAGs that were specifically enriched in L, such as L8_Streptomyces and L11_Myceligenerans, which displayed markedly higher relative abundances than others in this cluster. Group D represented MAGs more abundant in LW than in L, although without significant differences from other treatments. Groups E and F showed less consistent patterns. Group E contained MAGs with significant differences between treatments but no clear overall trend. Several members of this group, however, were abundant in C, P, and L, including halophilic taxa such as P10_Halomarina, C11_Actinopolyspora, P1_Haladaptatus, and C12_YR4-1 (an unclassified Balneolaceae). Group F comprised MAGs with no significant differences across treatments. Additionally, two MAGs assigned to taxon 70-9 (an unclassified Actinomycetota) were highly enriched in one replicate of PW (Fig. 5A).

Functional annotation using KEGG revealed that many MAGs shared broadly similar profiles at the KO identifier level, likely reflecting the large proportion of genes involved in conserved cellular processes (Fig. 5B). Nonetheless, clustering based on Jaccard distances explained 12.01% of the variance and showed significant differences, with MAGs in groups A and B displaying distinct metabolic signatures (PERMANOVA; p < 0.05, Fig. 5B). Taxonomic clustering further supported these patterns. For instance, halophilic archaea (P1, P10) and ammonia-oxidizing archaea (CW3, CW2) were positioned adjacently, while Actinobacteriota (e.g., order UBA5794) and unclassified Gammaproteobacteria (e.g., order UBA6522) were grouped as outliers (Fig. 5B).

Selection of lignocellulose-transforming microbial species
On average, 2,266 (± 875) genes were identified per MAG. Functional annotation using eggNOG-mapper revealed that the most abundant COG category was unknown function (Suppl. Data 4). In parallel, proteins from each MAG were annotated against the CAZy database (Suppl. Data 5), which resulted in 3,230 genes assigned to 190 families across the 74 MAGs. Analysis of CAZy families involved in lignocellulose degradation (Suppl. Data 5; Fig. 5C) identified 886 genes within the 74 MAGs, with most overrepresented in groups A and C (Fig. 5AC). Particularly, L8_Streptomyces, L11_Myceligenerans, and L14_Gimesia carried diverse repertoires of cellulases (e.g., AA10, GH16, GH5, GH6, GH48), xylanases (e.g., GH10, GH11, GH3, GH43, GH51), and esterases (e.g., CE1, CE4, CE7, CE14) (Suppl. Data 5). Three cellulose-degrading families (GH6, CBM4, and CBM16) were significantly enriched in L8_Streptomyces (ORA with Bonferroni correction; p < 0.05), while four families involved in cellulose and xylan deconstruction (AA10, GH11, GH42, GH48) were significantly enriched in L11_Myceligenerans (p < 0.05). Additionally, several MAGs in group A, such as LW3_Balneolaceae and LW13_Latescibacterales, contained a notable proportion of glycosyl hydrolases targeting plant polymers (e.g., GH10, GH43, GH51, GH67, GH53) (Suppl. Data 5).

Beyond polysaccharide degradation, we evaluated the presence and enrichment of 60 KEGG orthologs (KOs) linked to the catabolism of lignin-derived compounds. Within group A, genes encoding glutathione peroxidase (K00432), superoxide dismutase (K04564), and catalase (K03781) were significantly enriched in LW12_Alkalibacter, L2_Cohnella laeviribosi, and L1_Tuberibacillus calidus, respectively (ORA-Bonferroni correction; p < 0.05, Suppl. Data 6). The capacity to metabolize a range of lignin-derived aromatic monomers was observed across multiple MAGs (Suppl. Data 6; Fig. 5C). For example, genes encoding enzymes for catechol catabolism were enriched in C11 and LW11 (catA), PW1 (catC), and C3 and C16 (dmpB). All of these MAGs were affiliated with the Actinobacteriota phylum (ORA-Bonferroni correction; p < 0.05, Suppl. Data 6).

Screening of putative PET-depolymerizing enzymes within genomes
A total of 26 proteins derived from the 74 MAGs showed matches to enzymes in the PAZy database (> 50% amino acid similarity) (Suppl. Data 7). Particularly, two proteins from L8_Streptomyces and L11_Myceligenerans each had 26 hits against PAZy-derived proteins, displaying high similarity to PET40 (64% similarity, E-value 4.44E−101) and Cut190 (69% similarity, E-value 1.17E−118), respectively. In addition, LW4_Weizmannia coagulans contained four putative esterases with high similarity to HG-4 (99%), Est8-89 (58%), Est18-23 (61%), and Ces39-5 (60%), with E-values ranging from 2.91E−108 to 4.20E−229.

Lignocellulose amendment acts as a resource for putative PET-depolymerizing enzymes
To broaden the assessment of enzymatic potential beyond the assembled MAGs, we constructed a complete gene catalog—including genes from both binned and unbinned contigs—across all microcosm treatments. In total, 4,250,384 unique proteins were detected. This catalog was screened against the PAZy database using stringent thresholds (E-value ≤ 1E−15, ≥ 50% identity, and ≥ 80% coverage), resulting in 20–120 unique protein matches across treatments (Fig. 6A). Particularly, the lignocellulose-amended treatment (L) yielded the highest number of unique PAZymes, most of which were associated with PLA and PET transformation. In this treatment, approximately 60 putative PETases were detected, primarily affiliated with the Pseudomonadota, Actinomycetota, Bacillota, and Halobacteriota (Archaea) phyla (Fig. 6A). Interestingly, PETases from Halobacteriota were more frequently detected in dry treatments (P, C, and L). Within the entire gene catalog, only twelve unique proteins met ultra-stringent criteria (E-value ≤ 1E−20, ≥ 70% identity, and ≥ 90% coverage). Strikingly, all twelve were found in the L treatment. Among them, five were located in binned contigs affiliated with Micromonospora and Isoptericola, while four were affiliated with Marinobacter and Microbulbifer species (Fig. 6B).

Prediction of properties of putative PET-depolymerizing enzymes
A sequence-based comparative analysis was performed using two putative PETases identified in MAGs (L8_Streptomyces and L11_Myceligenerans), eleven putative PETases identified from the complete gene catalog (restricted to Pseudomonadota and Actinomycetota), and thirteen well-characterized PET-active enzymes (Suppl. Data 8). The SSN indicated that all enzymes clustered together, reflecting high alignment scores, although GlacPETase and HaloPETase1 were the most divergent (Fig. 6C). In the phylogenetic tree, these two enzymes grouped with Mic1 and Mic2 (Microbulbifer-derived putative PETases) (Fig. 6D). Interestingly, Mar1 and Mar2 proteins (from Marinobacter) clustered with HaloPETase5 and H. bauzanensis PETase, sharing ~72% amino acid similarity. The phylogenetic placement of Mic1/2 and Mar1/2 in the sequence landscape of PET-active enzymes from PAZy is shown in Suppl. Fig. 1. Among the 13 putative PETases, five lacked the functional M5 motif, including two from Micromonospora and two from Microbulbifer (Suppl. Data 9; Fig. 6D).

Clustering of predicted 3D structures showed that all type I terrestrial actinobacterial PETases grouped together, except two Micromonospora-derived enzymes (Fig. 7A). Additionally, both putative Microbulbifer-derived enzymes (Mic1 and Mic2) and HaloPETase1 displayed distinct structural features compared to Pseudomonadota homologs (Fig. 7A). Moreover, AI-based prediction of melting temperature (Tm) values via ProtScout revealed that three actinobacterial proteins had Tm values between 60 and 67 °C (Fig. 6E; Suppl. Data 10). Furthermore, ProtScout suggested that most putative PETases have an optimal activity at mesophilic conditions, >4% salinity, and near-neutral pH (5–9), except Mar2 from Marinobacter, which was predicted to function optimally under alkaline pH (9–14) (Suppl. Data 10; Fig. 7B). Multiple sequence and structure alignments confirmed the presence of catalytic triads (SDH), active sites, and disulfide bonds in all proteins (Suppl. Data 11; Fig 7B). For example, the catalytic triad of the L8_Streptomyces PETase was located at S138, D184, and H216, with an oxyanion hole at Y10 and M139, and a single disulfide bond between C249 and C265 (Fig. 6F). In contrast, both Marinobacter-derived proteins and PETase 403 (Ketobacter sp.) each contained three disulfide bonds (Suppl. Data 11; Fig. 7B).

Therefore, ESMfold-predicted 3D structures of Mic1 and Mic2 (Fig. 8A) as well as Mar1 and Mar2 (Fig. 8B) were aligned with crystal structures of their closest homologs: HaloPETase1 (PDB: 9hl5) for Mic1 and Mic2, and IsPETase (PDB: 5xjh) for Mar1 and Mar2 (see Fig. 7 for identification of related structural homologs). Overall, the active site architectures were qualitatively similar across all structures. However, root-mean-square deviation (RMSD) values of aligned backbone Cα atoms indicated greater structural divergence for Mic1 (1.50 Å) and Mic2 (1.23 Å) relative to HaloPETase1, compared with Mar1 (0.68 Å) and Mar2 (0.71 Å) relative to IsPETase (Fig. 8C and 8D). These global structural differences were also reflected locally at the predicted subsite I (classification according to [27]). While IsPETase homologs Mar1 and Mar2 showed highly similar residue occupations and spatial positioning of corresponding sidechains (Fig. 8D), HaloPETase1 homologs Mic1 and Mic2 displayed greater variation. Specifically, instead of T88, Q157, and Y180 in HaloPETase1, Mic1 and Mic2 contained Y, M, and a non-canonical H at the equivalent positions. At position L184, both candidates showed a conserved L in a spatially similar orientation. Notably, the positioning of sidechains H and L (corresponding to Y180 and L184 in HaloPETase1) in Mic1 and Mic2 may have been influenced by prediction uncertainty, as indicated by pLDDT scores below 80 (Fig. 8A). However, these lower pLDDT scores may reflect intrinsic conformational flexibility rather than prediction error alone.

Methods
Soil microcosm experiments
A laboratory-based microcosm experiment was conducted to induce shifts in the structure of mangrove soil microbial communities using inputs of lignocellulose (rice husk) and polyethylene terephthalate (PET) particles under two conditions: with and without seawater (collected near the soil sampling site). Before setting up the microcosms, rice husk was knife-milled through a 1-mm screen. PET particles (<5 mm; Eastar™ Copolyester GN071) and rice husk particles were washed twice with distilled water and 70% (vol/vol) ethanol, then air-dried at room temperature. A total of 25 soil samples (~500 g each) were collected from the top 10 cm of soil at the Barú peninsula (Cartagena de Indias, Colombia; 10°15′48.5′′N, 75°35′29.8′′W) (Fig. 1). To capture environmental heterogeneity, samples were taken at random intervals along a ~50 m transect extending inland from the shoreline. Soil and seawater sampling was conducted under permit number 2021123746-1-000 issued by the Autoridad Nacional de Licencias Ambientales (ANLA). Fresh mangrove soil samples were homogenized, sieved through a 2-mm mesh, and combined into a single bulk sample. Each experimental unit was prepared in a 100-mL Erlenmeyer flask containing 15 g of homogenized mangrove soil. An external carbon source (1.5 g), either rice husk (L treatment) or PET particles (P treatment), was then added at a final concentration of 10% w/w. Additionally, 15 mL of seawater was added to half of the experimental units, yielding the LW and PW treatments. Control microcosms with seawater (CW) or without seawater (C) received no additional carbon source. All microcosm treatments (L, LW, P, PW, C, and CW) were incubated at 30 °C for 30 and 90 days. Each treatment was replicated four times, resulting in a total of 48 experimental units (n = 48) (Fig. 1). After incubation, soils from the microcosms were homogenized and collected for DNA extraction using the DNeasy PowerSoil Kit (Qiagen, Maryland, USA), following the manufacturer’s instructions. DNA from the original (in-situ) homogenized soil was also extracted for comparison.

Microbial diversity analyses based on 16S rRNA gene and ITS2 amplicon sequencing
DNA extracted from soils was sent to Macrogen, Inc. (Seoul, Korea) for bacterial 16S rRNA gene amplicon sequencing of the V3–V4 region (primers Bakt-341F and Bakt-805R; [84]) and fungal internal transcribed spacer (ITS) 2 region (primers ITS3 and ITS4) using the Illumina MiSeq platform (300 bp paired-end reads). Raw sequences were processed with the DADA2 pipeline v1.26 (85). Briefly, primers were removed using Cutadapt (86), after which sequences were filtered and trimmed, errors were identified and corrected, and paired sequences were merged and dereplicated before taxonomy was assigned using the SILVA rRNA gene database (v.138) for bacteria and the UNITE database (v.10.5.2021) for fungi. The processed sequences were imported into PhyloSeq (v1.48.0), and non-fungal and non-bacterial reads were removed. Amplicon sequencing variant (ASV) tables were rarefied to 10,000 (16S rRNA gene) or 5,000 (ITS2) reads for downstream analyses, including comparisons of alpha and beta diversity metrics and statistical tests (e.g., PERMANOVA). All alpha diversity tests were conducted using ANOVA followed by Tukey’s HSD for pairwise comparisons (alpha = 0.05). ASVs associated with distinct treatment groups were then identified using the Boruta algorithm (87) and DESeq2 v1.5.0.2 (88). In particular, the Boruta algorithm was implemented in R and used to detect the features (i.e., ASVs) that explained variations between the treatment groups (i.e., seawater addition or lignocellulose addition). Briefly, bacterial ASVs indicative of treatment groups were modeled independently using random forest analysis with 1,000 trees followed by feature selection using the Boruta algorithm. The abundance of significant ASVs (e.g., importance > 2) were then z-score transformed and visualized with heat maps built using ggplot2 as is reported in (24). Differential abundance analysis was performed parametrically using the Wald test. The statistical significance of differentially abundant taxa was determined at p = 0.05 after FDR multiple comparisons corrections.

Additionally, nonrandom co-occurrence analyses were performed using SparCC (89). To this end, we selected bacterial and fungal ASVs with at least 30 sequences, which together accounted for more than 90% of the total reads. For each co-occurrence network, p-values were estimated from 99 permutations of randomly shuffled data tables. Only SparCC correlations with an absolute magnitude >0.7 or < −0.7 and statistical significance (p < 0.01) were retained for network construction. In the resulting networks, nodes represent ASVs, and edges denote strong and significant correlations. Network topology was assessed using metrics such as the number of nodes and edges, average degree, clustering coefficient, modularity, and network density. Network visualizations were generated in Gephi v0.1 (90).

Metagenome sequencing and read-based taxonomic classification
Three biological replicates from each treatment (i.e., L, LW, P, PW, C, and CW) incubated for 90 days (n = 18) were selected for metagenomic sequencing and analysis (Fig. 1). High-quality microbial DNA was sequenced at Novogene (Cambridge, UK) on the Illumina NovaSeq platform (2 × 150 bp; Illumina Inc., San Diego, CA, USA). The quality of raw reads was first assessed using FastQC v0.11.9 (91). Reads were then quality-filtered and trimmed with Trimmomatic v0.36 (minimum PHRED quality score of 30 and minimum length of 50 bases; [92]). Taxonomic profiling was performed on unassembled high-quality reads using Kraken2 v1.4.0 (93), with Bracken v2.9 (94) used to estimate relative abundances. The PlusPFP database (version 2024-09-04), which includes RefSeq entries for Archaea, Bacteria, Eukarya, viruses, and plasmids, served as the reference. Alpha and beta diversity metrics, including the Shannon index and Bray–Curtis dissimilarity, were calculated from the unassembled read-based taxonomic profiles using the vegan package in R (95). Complementary statistical analyses and visualizations were carried out in Python using the scipy (96) and scikit-bio v0.5.9 (97) libraries.

Reconstruction and analysis of metagenome-assembled genomes
For each treatment, high-quality and clean reads were co-assembled into contigs using the de novo assembler MEGAHIT v1.1.3 with the default k-mer range (98). Assembly quality was evaluated with QUAST v5.0.2 (99). Clean reads were mapped to contigs using Bowtie2 v2.4.5 (100), and contigs were binned with MaxBin2 v2.2.7 (101), MetaBAT2 v2.15 (102), and VAMB v4.1.3 (103), using default parameters and excluding contigs shorter than 1,000 bases. Bin refinement across treatments was performed using DAS Tool v1.1.4 (104), and bin quality was assessed with CheckM v1.2.1 (105). MAGs were then selected by retaining only medium- and high-quality bins (completeness ≥ 50% and contamination ≤ 10%), following the MIMAG standard (106). MAGs were taxonomically classified using GTDB-Tk v1.4.0 (107) using the GTDB R214 release (2023-04-28). Relative abundance of MAGs across samples was estimated with CoverM v0.7.0 (https://github.com/wwood/CoverM) by mapping quality-filtered reads against the entire MAG catalog. Since MAGs were recovered across all treatments rather than being exclusive to specific microcosms, a new grouping was established using co-abundance clustering. Kruskal–Wallis tests were first applied to relative abundance data to identify MAGs showing significant differences among treatments. Bonferroni correction was applied when conducting multiple comparisons. When significant effects were found, Dunn’s post hoc test was performed to assess pairwise differences. Finally, hierarchical clustering of MAGs was conducted using Ward’s linkage method to minimize intra-cluster variance. The resulting dendrogram defined clusters of co-abundant MAGs, which were subsequently used as categorical variables in downstream MAG-level analyses.

Functional annotation of genes derived from metagenome-assembled genomes
For each MAG, protein-coding sequences (CDSs) were predicted using the DFAST pipeline v1.2.6 (108), with Prodigal (109) employed for gene calling. Functional gene annotation was performed with eggNOG-mapper v2 using both the Clusters of Orthologous Groups (COG) and KEGG databases. To specifically evaluate lignocellulose-degrading potential, genes involved in lignin-transforming metabolic pathways were identified from KEGG based on curated gene sets reported by Díaz-García et al. (110). The presence of carbohydrate-active enzymes (CAZymes) was assessed using the dbCAN3 meta server (111). To assess the potential for plastic degradation, predicted CDSs were aligned against the PAZy database (112) using DIAMOND v0.9 (113) with a minimum amino acid identity threshold of 50%. To account for genome incompleteness, gene counts from COG, KEGG, CAZy, and PAZy annotations were normalized by dividing each value by the corresponding MAG completeness proportion, following the recommendations of Eisenhofer et al. (114). Functional divergence among MAGs was evaluated via principal coordinates analysis (PCoA) based on Jaccard distances of KEGG-annotated gene profiles (KO IDs), allowing visualization of functional clustering across MAG categories. PERMANOVA was performed to statistically test the significance of these functional differences. Results were visualized as ordination plots, with each point representing an individual MAG. Additionally, over-representation analysis (ORA) was conducted for genes associated with (hemi)cellulose and lignin transformations. Enrichment significance was determined using hypergeometric distribution tests with Bonferroni correction for multiple comparisons. All analyses and visualizations were carried out in Python, employing the scipy (96) and scikit-bio v0.5.9 (97) libraries.

Generation of a gene catalog from metagenomes and screening of plastic-active enzymes
To generate a comprehensive protein catalog from the microcosm-derived metagenomes, the open-source software VEBA was used (115). Briefly, raw sequences were co-assembled and binned using default workflows and parameters. Binned and unbinned contigs were taxonomically classified using GTDB-Tk (107) via VEBA pipeline. Gene prediction was performed with Prodigal (109), and the resulting gene models (GFF3 format) were modified to include gene and contig identifiers for downstream analyses. Genes were identified on both binned and unbinned contigs to ensure no information was lost, resulting in a complete gene catalog from all metagenomes. All predicted proteins were aligned against the PAZy database, which contains enzymes with experimentally validated activity on plastic polymers (112), using DIAMOND v0.9 (113). Two thresholds of amino acid sequence identity were applied: (i) high confidence (E-value ≤ 1E-15, ≥ 50% identity, and ≥ 80% coverage) and (ii) ultra-high confidence (E-value ≤ 1E-20, ≥ 70% identity, and ≥ 90% coverage). Unique protein counts were assigned to specific plastic polymers and normalized based on the total number of proteins in each treatment (L, LW, P, PW, C, and CW). Taxonomic affiliation of predicted plastic-active enzymes (PAZymes) located on binned contigs was determined using GTDB-Tk v1.4.0, while those on unbinned contigs were assigned based on the best NCBI BLAST hit.

Analysis of predictive PET-depolymerizing enzymes
Predicted PETases from the Pseudomonadota and Actinomycetota phyla that exhibited ultra-high similarity to enzymes in the PAZy database were analyzed in comparison with a set of known wild-type PET-active enzymes (e.g., PHL7, LCC, and IsPETase), including HaloPETase 1 and 5 (28) (Suppl. File 1), as well as two putative PETases detected in MAGs L8 and L11. Sequence Similarity Networks (SSNs) were generated using The Enzyme Function Initiative–Enzyme Similarity Tool (116) with an E-value threshold of 1E−10 and a sequence identity cutoff of 35% (nodes representing sequences sharing >35% identity were connected by edges, defining SSN clusters). The resulting SSN was visualized using Cytoscape (117). Multiple sequence alignment (MSA) was performed with MUSCLE (118), and phylogenetic trees were inferred using FastTree 2 (119) and visualized with iTOL (120). BLASTp was employed to compare sequences against the PETase from I. sakaiensis to identify the large functional M5 motif (73). Signal peptides were predicted using SignalP 6.0 (121) and removed prior to 3D structure prediction. For cases with no SignalP prediction, PrediSI (122) was used as an alternative. Three-dimensional structures were modeled via ESMFold, and corresponding sequence embeddings were extracted from the ESM-2 model (123). LE–His tags and predicted signal peptides were excised before docking. Molecular docking of a PET dimer substrate was performed using DiffDock-L (124) to generate plausible ligand–protein binding poses. Catalytic residues were identified by mapping conserved active-site positions from IsPETase onto the MSAs, guided by published structural data (27). Enzyme biophysical properties were ranked using the ProtScout Python package (https://github.com/new-atlantis-labs/ProtScout), incorporating AI-based models to predict melting temperature (Tm) via TemBERTure (125), optimal temperature, pH, and salinity via GeoPoc (126), and solubility via GATSol (127). Kinetic parameters (kcat and KM) for PET dimer hydrolysis were predicted using CatPred (128), a machine-learning tool trained on enzyme kinetics datasets. Classical physicochemical properties (e.g., molecular weight, isoelectric point) were computed using Biopython (129). Pairwise structural similarity was quantified with TM-align (130) to obtain TM-scores and RMSDs, which were then used as a distance metric for hierarchical clustering of the structures.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an excellent and timely manuscript that provides valuable insight into the emergence of halophilic plastic-degrading\nenzymes and the ecological processes driving microbial adaptation under saline and desiccating conditions. The\nexperimental design is solid and ecologically relevant, combining controlled microcosms with lignocellulose and\npolyethylene exposure to test selective pressures that mimic natural environments. The discovery of novel, potentially\nthermostable PETase-like enzymes highlights the metabolic versatility of extremophilic\ncommunities and will be of broad significance to environmental microbiology and biodegradation research.\nThe study is methodologically solid, with deep metagenomic sequencing (~65 million reads per sample; ~23.4 Gbp per\ntreatment) enabling both genome-resolved and unbinned functional analyses. The use of the VEBA pipeline and the\nPlusPFP database is a major strength, as it supports standardized and comprehensive protein catalog generation across\ndomains. However, the link between function and taxonomy remains unclear and, in its current form, not fully convincing. It\nappears that PETase-like sequences were identified from the unbinned protein catalog rather than\nwithin specific MAGs. If\nso, I would encourage the authors to clarify whether they explicitly searched for PETase genes within their MAGs and to\ndiscuss why these enzymes might not have been detected there (e.g., due to short contigs, plasmid localization, or\nlimitations of binning). This clarification is essential to assess whether the absence of PETases in MAGs reflects biological\nreality or methodological constraints. Additionally, a phylogenetic analysis of the PETase-like\nsequences, placed alongside\ntaxonomically validated reference PETases, would help verify their evolutionary relationships and possible genomic origins.\nSuch an analysis could substantially strengthen the functional–taxonomic linkage and support the claim that these enzymes\nrepresent novel or halophilic variants.\nThe integration of 16S rRNA amplicon and shotgun metagenomic datasets is another strength, though the analytical\nrelationship between them could be clarified. It would be helpful to specify whether 16S sequencing was used exclusively\nfor community structure and diversity analyses, or if it also informed metagenomic taxonomic profiles. Additionally, a brief\nexplanation of the Boruta algorithm (a random forest–based feature selection method) would improve methodological\ntransparency.\nAnother valuable but unexplored aspect is the eukaryotic fraction. Fungal sequences appear to follow similar ecological\npatterns to bacterial degraders, particularly under lignocellulose exposure and reduced salinity. Given the well-documented\nrole of fungi in producing cutinases, esterases, and oxidoreductases, a short discussion acknowledging this potential\ncontribution would broaden the ecological scope of the study.\nOverall, this is a high-quality and well-executed piece of research. The conclusions are supported by the data, and the\nmethodology meets current standards in environmental metagenomics. The main areas for improvement are (i) a clearer\nexplanation of how VEBA, PlusPFP, and MAG reconstruction were used in tandem; (ii) explicit acknowledgment of the\nlimitations of the functional–taxonomic linkage; and (iii) a brief perspective on the fungal component. Addressing these\npoints would further strengthen an already impressive and impactful contribution.", "reviewer_2": "Comments to Pena-Valencia et al.\nThe manuscript entitled “Lignocellulose-mediated selection of halophilic PET-degrading enzymes from mangrove soil” by\nPena-Valencia et al. reported prokaryotic and fungal community perturbations in microcosm incubations induced by the\naddition of seawater, PET, and lignin. The authors used marker-gene amplicon and metagenome sequencing to\ncharacterize the community changes (particularly PET enzyme-containing species enrichment) in different treatments. I am\nnot a PET expert and cannot judge the parts related to PET enzymes and structures. I don’t know how big biotechnological\napplication potentials these findings have. I apologize first if my criticisms below delay the publication of this manuscript.\nI find it lack focal points from a microbiologist’s perspective. It covers many aspects, but lack deep insights on most of them.\nSequencing data alone cannot decipher the complex communities and their responses to perturbations in mangrove soils. I\nhave the following two suggestions for the authors to consider collecting some wet lab data to reinforce their findings.\nIt is clear that some lineages became more abundant than others, which can be due to multiple reasons. Is it possible that\nthey are just more resistant to the imposed unfavorable conditions relative to others? I can see some of the discussed\norganisms seems to be quite cosmopolitan and can be seen in many human-affected environments. Their abundance\nchanges over time could be the key to really reveal the population dynamics (growth/decay) of each discussed organismal\ngroup.\nA lot of PET enzyme findings were presented and discussed, based on the dry lab sequencing data. Are their any data to\nsupport that those organisms are really responding to the imposed PET/lignin additions in the incubations, like transcripts?\nMoreover, is it possible to measure the concentrations of the added PET/lignin or their degradation products over time to\nprove the activities of the bulk communities? Finally, is it possible to have some PET-enzyme assays to check if they have\nthe assumed activities or not."}, "author_response": "Responses to the reviewer’s comments are as described below (in blue):\nReviewer #1:\nThis is an excellent and timely manuscript that provides valuable insight into the emergence of\nhalophilic plastic-degrading enzymes and the ecological processes driving microbial adaptation\nunder saline and desiccating conditions. The experimental design is solid and ecologically\nrelevant, combining controlled microcosms with lignocellulose and polyethylene exposure to test\nselective pressures that mimic natural environments. The discovery of novel, potentially\nthermostable PETase-like enzymes highlights the metabolic versatility of extremophilic\ncommunities and will be of broad significance to environmental microbiology and biodegradation\nresearch.\nR/ Thank you for the positive assessment of our article. We are glad that you have found it solid,\ntimely, and interesting.\nThe study is methodologically solid, with deep metagenomic sequencing (~65 million reads per\nsample; ~23.4 Gbp per treatment) enabling both genome-resolved and unbinned functional\nanalyses. The use of the VEBA pipeline and the PlusPFP database is a major strength, as it supports\nstandardized and comprehensive protein catalog generation across domains. However, the link\nbetween function and taxonomy remains unclear and, in its current form, not fully convincing. It\nappears that PETase-like sequences were identified from the unbinned protein catalog rather than\nwithin specific MAGs. If so, I would encourage the authors to clarify whether they explicitly\nsearched for PETase genes within their MAGs and to discuss why these enzymes might not have\nbeen detected there (e.g., due to short contigs, plasmid localization, or limitations of binning). This\nclarification is essential to assess whether the absence of PETases in MAGs reflects biological\nreality or methodological constraints. Additionally, a phylogenetic analysis of the PETase-like\nsequences, placed alongside taxonomically validated reference PETases, would help verify their\nevolutionary relationships and possible genomic origins. Such an analysis could substantially\nstrengthen the functional–taxonomic linkage and support the claim that these enzymes represent\nnovel or halophilic variants.\nR/ Indeed, the metagenomic information obtained from MAGs can be constrained by the binning\nprocess, where multiple low-quality bins, unbinned contigs, and short sequences are excluded in\nsubsequent analyses (e.g., functional gene annotation). This becomes an important issue to be\nconsidered in environmental samples containing high levels of microbial diversity, which cause a\nlow recovery of high-quality MAGs (see discussion in lines 378-391). Thus, to circumvent this\nissue, our search for predictive PET hydrolases (PETases) was performed in high-quality MAGs\nand in the full gene catalog obtained -via the VEBA pipeline- from binned and unbinned contigs\n(see lines 262-266). As we mentioned in the text, MAGs provide a reasonable taxonomic\nresolution, sufficient to establish an accurate link between taxa and the detected enzymes (see lines\n382-384). In our study, the majority of twelve putative PETases-like enzymes (except L8 and L11)\nwere detected within unbinned contigs or low-quality bins and not within the set of high-quality\nMAGs. This suggests that some of these enzymes are likely within species present at low\nabundance in our samples, which limits their MAG recovery. In our study, this issue was surpassed\nby using both gene and genome-resolved metagenomics. We clarified this point in the revised text\n(see lines 406-411). In addition, taxonomy affiliation of binned and unbinned contigs was also\ncarried out via the VEBA pipeline (information included in methods, see lines 601-602), this helps\nto clarify and support functional–taxonomic linkage.\nMoreover, comprehensive phylogenetic analyses were performed in our study. For instance, in\nFig. 6D, we built a phylogenetic tree using all putative PETases detected in MAGs and the full\ngene catalog from treatment L, in addition to thirteen prototypical and validated PET-active\nreference enzymes (such as LCC, IsPETase, and PHL7) deposited in the PAZy database\n(https://onlinelibrary.wiley.com/doi/10.1002/prot.26325). To complement this analysis, the same\nset of enzymes was used to build a phylogenetic tree using 3D structural features. Additionally, a\ncomplete phylogenetic tree was reconstructed using all PET-active enzymes available in the PAZy\ndatabase alongside four putative PETases from Marinobacter (Mar1 and Mar2) and Microbulbifier\n(Mic1 and Mic2) (see Suppl. Fig. 2, lines 286-287). All these computational analyses allowed us\nto gain insights into the evolutionary relationship, novelty, and taxonomic origin of predictive\nPETases.\nmes lacking the functional M5 motif\nThe integration of 16S rRNA amplicon and shotgun metagenomic datasets is another strength,\nthough the analytical relationship between them could be clarified. It would be helpful to specify\nwhether 16S sequencing was used exclusively for community structure and diversity analyses, or\nif it also informed metagenomic taxonomic profiles. Additionally, a brief explanation of the Boruta\nalgorithm (a random forest–based feature selection method) would improve methodological\ntransparency.\nR/ For clarity, whole metagenomic sequencing was performed only at the end of the microcosm\nexperiment (i.e., after 90 days of incubation). This decision was based on the information obtained\nvia bacterial 16S rRNA gene amplicon sequencing, where drastic changes in the initial microbial\ncommunities' structure were observed (see Fig. 2). These changes were mostly associated with\nselection caused by desiccation and salinity. This information is now provided in the revised text\n(see lines 163-166).\nThe Boruta algorithm is a machine learning algorithm used for feature selection in highly\ndimensional datasets, including microbiome data (see\nhttps://pmc.ncbi.nlm.nih.gov/articles/PMC6433899/). The approach allows for the identification\nof features that account for differences in microbial community composition or functional profiles\nacross treatment groups. Briefly, the Boruta algorithm identifies important features (i.e., ASVs)\nby comparing the importance of real predictors to that of randomized “shadow” features. These\nshadow features are generated by permuting the original data across samples. In doing so, this\npermutation removes the real structure and the associations of treatments with response variables\n(i.e., ASVs). Importance values of the permuted data are generated through Random Forest and\nthen compared to the real data. A feature is only considered important if its Random Forest\nimportance score is consistently greater than that of the best-performing randomized variable. In\nother words, for an ASV to be important, it must do a better job of predicting the treatment group\nthan the best-performing randomized value. This approach ensures that the selected features are\nmore informative than what could be expected by random chance. Shadow feature comparison is\nwell-suited for microbiome data (i.e., high-dimensional and sparse), where many ASVs may be\nweakly predictive due to random chance as opposed to real treatment associations.\nThe Boruta algorithm was implemented in R (see:\nhttps://www.jstatsoft.org/article/view/v036i11/0) and used to detect the features (i.e., ASVs) that\nexplained variations between the treatment groups (i.e., seawater addition or lignocellulose\naddition). Briefly, bacterial ASVs indicative of treatment groups were modeled independently\nusing random forest analysis with 1,000 trees, followed by feature selection using the Boruta\nalgorithm. The abundance of significant ASVs (e.g., importance > 2) was then z-score transformed\nand visualized with heat maps built using ggplot2 (as is reported in\nhttps://www.sciencedirect.com/science/article/pii/S0167779924002427#bbb0570). This\ninformation is now included in the methods section (see lines 521-527).\nAnother valuable but unexplored aspect is the eukaryotic fraction. Fungal sequences appear to\nfollow similar ecological patterns to bacterial degraders, particularly under lignocellulose exposure\nand reduced salinity. Given the well-documented role of fungi in producing cutinases, esterases,\nand oxidoreductases, a short discussion acknowledging this potential contribution would broaden\nthe ecological scope of the study.\nR/ We agree that fungal species can harbor potentially multiple biotechnologically interesting\nenzymes. However, in our study, bacterial-derived sequences dominated in the metagenome data\n(~ 95% in most of the samples) with a very low proportion of unassembled sequences (~ 0.2%)\nassigned to fungi via Kraken2. We did not consider exploring these low-abundant sequences and\nrevised this information in the main text (see lines 391-395).\nOverall, this is a high-quality and well-executed piece of research. The conclusions are supported\nby the data, and the methodology meets current standards in environmental metagenomics. The\nmain areas for improvement are (i) a clearer explanation of how VEBA, PlusPFP, and MAG\nreconstruction were used in tandem; (ii) explicit acknowledgment of the limitations of the\nfunctional–taxonomic linkage; and (iii) a brief perspective on the fungal component. Addressing\nthese points would further strengthen an already impressive and impactful contribution.\nR/ Thanks for your comments. First, Kraken/Bracken-PlusPFP was used to taxonomically classify\nunassembled sequences, the VEBA pipeline was used to develop a full gene catalog from binned\nand unbinned contigs, and MAG reconstruction was carried out to obtain taxonomic and functional\ninsights into the most abundant microorganisms within the soil microcosms. Second, the\ntaxonomic origin of putative PETases was determined by the taxonomic affiliation of binned or\nunbinned contigs (via VEBA) where genes were detected. This taxonomic assignment was\nsupported by amino acid sequence analyses and 3D structure comparison along prototypical PET-\nactive enzymes. Third, the fungal component and its potential were mentioned in the revised text.\nReviewer #2:\nComments to Pena-Valencia et al.\nThe manuscript entitled “Lignocellulose-mediated selection of halophilic PET-degrading enzymes\nfrom mangrove soil” by Pena-Valencia et al. reported prokaryotic and fungal community\nperturbations in microcosm incubations induced by the addition of seawater, PET, and lignin. The\nauthors used marker-gene amplicon and metagenome sequencing to characterize the community\nchanges (particularly PET enzyme-containing species enrichment) in different treatments. I am not\na PET expert and cannot judge the parts related to PET enzymes and structures. I don’t know how\nbig biotechnological application potentials these findings have. I apologize first if my criticisms\nbelow delay the publication of this manuscript.\nR/ Thank you for all the comments provided. For clarity, our study is located between the microbial\necology-enzyme bioprospecting/discovery zone. It provides (i) diversity- and genomic-based\ninsights that advance our understanding of mangrove soil microbiomes under changing\nenvironmental conditions, and (ii) a conceptual and methodological framework to improve the\nselection of putative PET-degrading enzymes from environmental microbiomes. Importantly, this\nstudy does not aim to demonstrate or validate the industrial applicability of the identified enzymes.\nRather, based on amino acid sequence analyses, predicted 3D structures, and AI-assisted\npredictions, we reported enzymatic features (e.g., halotolerance and thermostability) that suggest\nthese candidate PETases may be promising targets for downstream experimental characterization,\nbenchmarking, and evaluation under industrially relevant conditions.\nI find it lack focal points from a microbiologist’s perspective. It covers many aspects, but lack\ndeep insights on most of them. Sequencing data alone cannot decipher the complex communities\nand their responses to perturbations in mangrove soils. I have the following two suggestions for\nthe authors to consider collecting some wet lab data to reinforce their findings.\nIt is clear that some lineages became more abundant than others, which can be due to multiple\nreasons. Is it possible that they are just more resistant to the imposed unfavorable conditions\nrelative to others? I can see some of the discussed organisms seems to be quite cosmopolitan and\ncan be seen in many human-affected environments. Their abundance changes over time could be\nthe key to really reveal the population dynamics (growth/decay) of each discussed organismal\ngroup.\nR/ Thank you for this positive feedback provided and to highlight potential issues under the view\nof microbiologists. As specifically discussed in our article, the absence of seawater, the inclusion\nof lignocellulose, and the often high salinity conditions caused by desiccation were the major\ndrivers that selectively impacted the microbial community structure in our microcosm system (see\nlines 341-354). These factors operating over time resulted in the selection of resistant\nmicroorganisms (e.g., spore-forming, halotolerant, and lignocellulolytic ones) after 90 days of\nincubation in treatment L, increasing their relative abundance compared with others. As mentioned\nby the reviewer, some of these selected taxa belong to cosmopolitan halophilic groups (e.g.,\nBacillacea), although we have also found specialized taxa thriving under these conditions in our\nmicrocosm system (e.g., Micromonosporaceae). Probably, extensive sampling over time can\nprovide information about population dynamics; however, this goal was not initially considered in\nour experimental design. Based on your comment, we acknowledged in the revised text that\nalternative techniques such as metatranscriptomics can provide insight in microbial activity and\nalso population dynamics after perturbation events (see lines 371-376).\nA lot of PET enzyme findings were presented and discussed, based on the dry lab sequencing data.\nAre their any data to support that those organisms are really responding to the imposed PET/lignin\nadditions in the incubations, like transcripts? Moreover, is it possible to measure the concentrations\nof the added PET/lignin or their degradation products over time to prove the activities of the bulk\ncommunities? Finally, is it possible to have some PET-enzyme assays to check if they have\nthe assumed activities or not.\nR/ We acknowledge the importance of the points raised by the reviewer. We consider that the\nmetagenomic-derived (dry-lab sequencing data) results are sufficiently robust and support the\nclaims and arguments developed in the main text.\nIn our study, predictive PETases were enriched in the metagenomes from lignocellulose-added\nmicrocosm (L treatment) compared to controls (Fig. 6A). Moreover, specific bacterial taxa (e.g.,\nMicromonospora, Streptomyces, and Microbulbifer) that hold these putative PETases were\nprimarily responding to the input of rice husk, but not to the input of PET. In this regard, these\ntaxa are known for their capacity to degrade plant biomass, and they can be metabolically active\nwithin the microcosms. Moreover, rice husk material was almost completely consumed (see photos\nin Fig. 1B), and the remaining plant biomass (after 90 days of incubation) would be the hardest\nfraction to be degraded (i.e., lignin). We have mentioned this in the revised text (see lines 336-\n338). Likewise, the amendment with PET particles did not result in detectable effects on the\nmicrobial communities compared to controls (Fig. 2B). In addition, depolymerization of PET (and\nrelease and consumption of PET-derived monomers) is unlikely to occur in the established\nmicrocosm conditions. As we mentioned in the text, PET-added microcosms did not favor the\nselection of PET-degrading taxa, likely due to the recalcitrance of PET to depolymerization, the\navailability of alternative carbon sources, or harsh environmental conditions that constrain PET\ndegradation (see lines 397-400).\nIn the revised text, we also acknowledge that metatranscriptomic data can be used to validate the\nmetabolic activity of the enriched taxa at the time of sampling (see lines 371-376). However, our\ndual-approach (i.e., gene-centric and genome-resolved metagenomics) was useful to explore the\nmetabolic potential of microbial taxa within complex community contexts. Moreover, as the\nmicrocosm experiment was carried out in 2021, metatranscriptome sequencing, analysis of\ntranscripts, and measurement of lignin content from the same flasks is now technically impossible.\nWe carefully acknowledge some of these limitations in our study, see lines 371-376.\nLast, the predictive enzymatic activity of the twelve putative PETases is strongly supported by the\npresence of catalytic features, including the catalytic triad, disulfide bonds, the M5 motif, as well\nas phylogenetic placement, amino acid sequence similarity, and predicted 3D structural similarity\nto previously characterized PET-active enzymes. Comprehensive wet-lab validation of enzymatic\nactivity for these twelve candidates PETases would involve complex experimental processes,\nwould be time-consuming, and can face multiple experimental bottlenecks (e.g., heterologous\nexpression, protein solubility issues, and enzymatic tests). In addition, it might fall outside the\nscope of the present study, as detailed enzymatic activity characterization, benchmarking, and\nassessments of industrial applicability were not initially considered in our main aims (see aims in\nlines 113-119), Nevertheless, we explicitly acknowledge this next step and highlight it as a future\nresearch direction in the revised text (see conclusion section, lines 473-476)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfactorily addressed my concerns.", "reviewer_2": "The authors had mostly addressed my questions and concerns."}, "author_response": ""}]

Decision letter