Metagenomic investigation of bacterial laccases in a straw-amended soil

Soil sampling and DNA extraction

Soil samples were collected from an experimental field located at the Agricultural Station of the Chinese Academy of Agricultural Sciences (39°36′N, 116°36′E), situated at an elevation of 18 m above sea level in Langfang, Hebei, China. The experiment site, soil taxonomy and soil properties were described in more detail in the previous study (Yu et al., 2018). The field follows a winter wheat-summer maize rotation cropping system. Calcium superphosphate (150 kg P₂O₅ ha⁻¹) and potassium sulphate (75 kg K₂O ha⁻¹) were applied before tillage. After each harvest, the total crop residue was directly shredded and incorporated into the soil through rotary tillage, to a depth of approximately 20 cm. The quantities of wheat (44.7% of total C, 0.76% total of N) and maize (38.6% of total C, 0.92% of total N) straw incorporated were approximately 3,000 kg ha⁻¹ and 8,000 kg ha⁻¹, respectively. This straw amendment practice began in 2010, with three replicate plots of 67 m²each.

Soil samples were collected in early November, 2016. In each replicate plot, five soil cores (3.0 cm in diameter) were randomly extracted from the plow layer and combined into a composite sample, resulting in a total of three soil samples. Each sample was sealed in a sterile plastic bag and transported to the laboratory on ice. After homogenization, the samples were sieved through a 2.0 mm mesh and stored at −80 °C until further analysis.

Total DNA was extracted from each soil sample using the PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA), following the manufacturer’s protocol with minor modifications (Yu et al., 2018). The extracted DNA was initially assessed using 1.0% agarose gel electrophoresis, and its purity was measured using a NanoDrop ND-2000 UV-VIS spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). In addition, its concentration was determined using the Qubit dsDNA Assay Kit and a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Finally, the DNAs extracted from the three soil samples were mixed at the same concentration in equal volumes for metagenomic sequencing.

Metagenomic sequencing and analysis

The pooled soil DNA was utilized for metagenomic library preparation and sequencing. Metagenomic libraries were constructed using the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). Briefly, the DNA was fragmented to a target size of 350 bp by sonication. The resulting fragments were end-repaired, ligated to full-length adaptors, and subsequently used to create paired-end sequencing libraries. High-throughput sequencing was performed on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) at Novogene (Beijing, China), generating approximately 14.11 Gb of raw data. Low-quality reads containing more than 40 nucleotides with a quality score below 38, or more than 10 ambiguous bases (N), were filtered out. Additionally, adapter sequences were removed. After filtering, a total of 14.09 Gb of clean data was retained. These data were deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA420606. The clean reads were assembled using MEGAHIT v1.1.1 with the “–presets meta-large” parameter, producing scaffolds (Li et al., 2015). These scaffolds were split at N-joins to generate scaftigs (continuous sequences without ambiguous bases). Only scaftigs ≥ 500 nt were retained for further open reading frame (ORF) prediction using MetaGeneMark v2.10 (Zhu, Lomsadze & Borodovsky, 2010). ORFs shorter than 100 bp were discarded, and the remaining sequences were clustered using CD-HIT v4.5.8 at a threshold of 95% identity and 90% coverage to construct a non-redundant gene catalog (Fu et al., 2012). The longest sequence within each cluster was selected as the representative for functional annotation. To assess gene abundance, the clean reads were realigned to the gene catalog using SoapAligner v2.21, with the parameters: -m 200, -x 400, identity ≥ 95% (Li et al., 2009). Genes with fewer than two mapped reads were excluded to minimize erroneous identifications. The relative abundance of each gene was calculated based on the number of reads mapped to the gene and the gene’s length.

For taxonomic assignment, predicted genes were aligned against the MicroNR dataset, which consists of microbial reference sequences from the NCBI non-redundant (NCBI-NR) database (Version: 2015-12-04), using DIAMOND v0.8.22 (BLASTp, e-value ≤ 1e−5) (Buchfink, Xie & Huson, 2015). Significant matches for each gene were defined as those with e-values ≤ 10 times the e-value of the top hit. Taxonomic classifications were determined using the lowest common ancestor algorithm implemented in MEGAN4 (Huson et al., 2011). The relative abundance of each taxonomic group was calculated by summing the relative abundances of the corresponding genes. For functional classification, DIAMOND v0.8.22 was used to align predicted genes with sequences in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Carbohydrate-Active Enzymes (CAZy) databases (BLASTp, e-value ≤ 1e−5). The best hits were selected for functional annotation.

Phylogenetic analysis of bacterial laccases

The deduced amino acid sequences from the predicted ORFs were screened to identify putative bacterial laccases using previously established profile hidden Markov models (pHMMs) (Ausec et al., 2011). Only sequences with an e-value below 1e−10 were retained. These sequences were manually validated against the NCBI Conserved Domain Database to confirm the presence of the characteristic cupredoxin domain (Dandare et al., 2019). The regions spanning between copper-binding regions I (cbr I) and II (cbr II) of the identified bacterial laccases were aligned with known sequences from the GenBank database using ClustalX2. A phylogenetic tree was constructed using the neighbor-joining method in MEGA v7.0, with bootstrap values calculated from 1,000 replications to assess the robustness of the tree topology (Kumar, Stecher & Tamura, 2016).

Cloning of bacterial laccase genes from soil DNA and sequence analysis

A subset of bacterial laccase genes was randomly selected for direct PCR amplification using soil DNA as the template. The amplified products were purified using the TIANgel Midi Purification Kit (Tiangen, Beijing, China) and subsequently cloned into pEASY-Blunt simple cloning vectors (TransGen, Beijing, China) for sequencing. The sequencing results were compared to the original sequences using BLAST for validation. Full-length cloning was performed on five selected laccase gene fragments. The 5′ and 3′ flanking regions were obtained using modified thermal asymmetric interlaced PCR TAIL-PCR (mTAIL-PCR), and these sequences were assembled with the known portions of the gene (Huang et al., 2010). The ORFs were identified using ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder), and sequence identities were verified using BLAST. Signal peptides were predicted using SignalP 5.0 (http://www.cbs.dtu.dk/services/SignalP/), and multiple sequence alignments were generated using ClustalX2.

Heterologous expression and purification of lacS1 in Escherichia coli (E. coli)

The ORF of lacS1, excluding the signal peptide-coding region, was amplified from the soil metagenome using primers specifically designed for this study: ccatgg ATGCAACGCCGATTAGAAAGCTGTCG (italicized Nco I site) and gcggccgc CGCAACGCGCACGAGCGC (italicized Not I site). The PCR product was purified and cloned into a pEASY-Blunt simple cloning vector for sequencing. After verification, the lacS1 fragment was subcloned into the pET-28a(+) expression vector (Novagen; Merck KGaA, Darmstadt, Germany). The recombinant plasmid pET-lacS1 was then transformed into E. coli BL21 (DE3) competent cells.

Transformed E. coli BL21 (DE3) cells harboring pET28a-lacS1 were cultured at 37 °C in 1 L of Luria-Bertani (LB) medium containing 50 µg/mL kanamycin with agitation at 200 rpm. When the culture reached an OD₆₀₀ of 0.5–0.6, enzyme expression was induced with 0.1 mM isopropyl-β-D-thiogalactoside and 0.25 mM CuSO₄. To maximize the yield of soluble and active protein, agitation was stopped after 4 h of induction (28 °C, 150 rpm), and the cells were incubated under microaerobic conditions for an additional 20 h. The cells were harvested by centrifugation at 11,000 rpm (4 °C) for 15 min. The resulting cell pellets were resuspended in cold 20 mM Tris–HCl buffer (pH 7.6) containing 500 mM NaCl and 10% glycerol. Cell disruption was performed by sonication on ice using an ultrasonic cell disruptor (Scientz, Ningbo, China) with 3-s pulses and 5-s intervals. After centrifugation at 11,000× g for 30 min, the supernatant was subjected to affinity chromatography using a Ni²⁺-nitrilotriacetic acid (NTA) chelating column (Qiagen, Hilden, Germany). Elution was carried out using a stepwise gradient of imidazole (0, 20, 40, 60, 80, 100, 200, 300, and 500 mM) in 20 mM Tris–HCl, 500 mM NaCl, and 10% glycerol (pH 7.6). Fractions exhibiting enzymatic activity were collected, and the purity and molecular mass of the enzyme were assessed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE).

Characterization of lacS1

Laccase activity was measured using 2,2′-Azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) (ABTS; Sigma-Aldrich, St. Louis, MO, USA) as the substrate. The reaction mixture contained 50 µL of appropriately diluted enzyme solution and 950 µL of 0.5 mM ABTS (ɛ₄₂₀ = 36,000 M⁻¹ cm⁻¹) in 50 mM sodium citrate buffer (pH 5.0) supplemented with 100 µM CuSO₄. The mixture was incubated at 50 °C for 5 min, followed by rapid cooling in an ice-water bath for 30 s to stop the reaction. Absorbance at 420 nm was measured spectrophotometrically. One unit (U) of laccase activity was defined as the amount of enzyme required to oxidize one µmol of ABTS per min. All assays were performed in triplicate.

The optimal pH for the purified recombinant lacS1 was determined by assaying laccase activity at 50 °C in 50 mM sodium citrate buffer (pH 3.0–6.0) and Na₂HPO₄-KH₂PO₄ buffer (pH 6.0–8.0). pH stability was assessed by measuring residual enzyme activity under standard conditions after incubation at various pH values (3.0–9.0) for 1 h at 37 °C. The buffers used for stability assays included 50 mM sodium citrate for pH 3.0–6.0, Na₂HPO₄-KH₂PO₄ for pH 6.0–8.0, and Tris–HCl for pH 8.0–9.0. Optimal temperature was determined by conducting the reaction at various temperatures (30–70 °C) in 50 mM sodium citrate buffer (pH 5.0). Thermal stability was evaluated by pre-incubating the enzyme at 40 °C and 50 °C for different durations, followed by activity measurements under standard conditions.

To assess the effect of CuSO₄ on lacS1 activity, the purified enzyme was incubated with ABTS at the optimal temperature for 5 min in 50 mM sodium citrate buffer (pH 5.0) in the presence of 0–0.5 mM CuSO₄. The effect of NaCl was similarly tested under standard conditions with varying concentrations.

Overview of the metagenomic data analysis

To uncover bacterial laccases in the straw-amended soil, metagenomic sequencing was performed using the HiSeq 2500 platform. After quality control, approximately 99.86% of the reads were used for assembly, generating a total of 416,112 scaftigs. MetaGeneMark predicted 329,007 non-redundant ORFs from these sequences (Table S1).

Taxonomic analysis revealed that 308,011 (93.62%) of the predicted ORFs could be annotated using the microNR database. Bacteria were the most dominant group (88.11%), followed by Archaea (4.29%), Eukaryotes (0.09%), and Viruses (0.05%). At the phylum level, the most abundant bacterial groups (relative abundance ≥ 1%) were Proteobacteria, Actinobacteria, Acidobacteria, Gemmatimonadetes, Bacteroidetes, Thaumarchaeota, Nitrospirae, Chloroflexi, Cyanobacteria, Firmicutes, Planctomycetes, and Verrucomicrobia (Fig. 1).

Functional analysis of the microbial communities in the straw-amended soil was performed by comparing the predicted proteins from the metagenomic dataset to the KEGG database. The KEGG annotations indicated that carbohydrate metabolism was the most abundant functional category (Fig. S1). Additionally, the predicted genes were analyzed against the CAZy database to gain insights into the microbial degradation of plant polymers in the soil. Most of the identified CAZymes were categorized as glycoside hydrolases (GHs, 37.8%), followed by glycosyltransferases (GTs, 33.9%), carbohydrate-binding modules (CBMs, 23.9%), carbohydrate esterases (CEs, 6.8%), auxiliary activities (AAs, 4.2%), and polysaccharide lyases (PLs, 0.8%) (Table 1). The CAZyme repertoire consisted of 4,727 candidate GHs from 92 families, 4,221 GTs from 43 families, 2,980 CBMs from 38 families, 853 CEs from 14 families, 522 AAs from nine families, and 106 PLs from 12 families. These results suggest that the microbial consortia in the straw-amended soil possess substantial capabilities for plant biomass degradation.

Diversity and phylogenetic distribution of bacterial laccases in the straw-amended soil

To identify putative bacterial laccase-coding genes in the straw-amended soil metagenome, pHMMs were employed to search the deduced amino acid sequences of all predicted ORFs (Ausec et al., 2011). This analysis identified 322 putative bacterial laccases, including both two- and three-domain enzymes (Table 2). Notably, 85.4% of these putative laccases were retrieved using multiple models, likely due to the incomplete nature of many sequences.

To validate the metagenomic sequencing and assembly, 14 putative bacterial laccase genes were randomly selected for PCR amplification and sequencing. Among these, 12 PCR products were obtained with the expected sizes. Further sequencing showed that only one fragment had no significant identity with the corresponding predicted gene, while the remaining 11 fragments exhibited sequence identities above 96%, except for two sequences with 86% and 92% identity. These results suggest that most of the predicted laccase genes represent authentic laccase genes in the straw-amended soil microorganisms.

The amino acid sequences of the identified laccases were aligned with known sequences in the NCBI-NR, CAZy, NCBI environmental (NCBI-ENV), and Swiss-Prot databases to assess their similarity with known proteins (Fig. 2). BLASTp analysis against the NCBI-NR database showed that the sequence identities of the 322 putative bacterial laccases ranged from 41% to 97%, with an average identity of 70.22%. Approximately 24.84% of these sequences were most similar to proteins annotated as “hypothetical proteins” in NCBI-NR. None of the sequences exhibited high similarity (>95% sequence identity) to any entry in the CAZy database, suggesting that none of these sequences had been previously deposited in CAZy (Hess et al., 2011; Zhou et al., 2017). Additionally, 45 sequences with less than 30% identity were considered novel (Zhou et al., 2017; Pearson, 2013). Furthermore, only 4.66% of the putative laccases shared more than 75% identity with sequences in the NCBI-ENV database, indicating that these enzymes have not been widely identified in previous metagenomic studies (Hess et al., 2011; Zhou et al., 2017). Lastly, 48 sequences exhibited less than 30% identity to any known proteins in Swiss-Prot, further indicating that their activity has not been biochemically verified (Zhou et al., 2017). Collectively, these results highlight the novelty of the bacterial laccases identified in the straw-amended soil.

The phylogenetic origins of the identified laccase genes were further explored through taxonomic classification (Fig. 3). These genes were found to be unevenly distributed across various phyla, with the majority originating from Proteobacteria (37.21%), Actinobacteria (16.32%), and Gemmatimonadetes (13.50%). Additional phyla contributing to laccase production included Thaumarchaeota (7.47%), Bacteroidetes (3.43%), Firmicutes (2.81%), Chloroflexi (1.53%), Crenarchaeota (1.08%), and Armatimonadetes (1.05%). Other phyla such as Acidobacteria, Cyanobacteria, Planctomycetes, Acetothermia, and Nitrospirae were also detected, albeit in lower relative abundances. At the order level, the genes were distributed across 42 different orders, with 17 orders representing over 1% of the total abundance (Fig. 3). The highest number of laccase-coding genes was affiliated with Gemmatimonadales (13.50%), followed by Geodermatophilales (6.86%), Nitrososphaerales (4.97%), and Burkholderiales (4.37%) (Fig. 3).

To further explore the diversity and relationship of bacterial laccases in the straw-amended soil, amino acid sequences spanning between cbr I and cbr II of the putative laccases were used to construct a phylogenetic tree. Therefore, 74 putative laccases were selected, and 29 representative sequences from 20 genera and environmental samples were used as references. The bacterial laccase fragments were separated into seven distinct clusters (A–G) (Fig. 4). The laccase fragments from the straw-amended soil were found in all seven clusters, underscoring their substantial diversity. Cluster A and Cluster E were the most populated, containing 23 and over one-third of the bacterial laccase fragments from the soil, respectively. Cluster B contained three laccases with no known references, suggesting high novelty. Clusters C and F had two bacterial laccase fragments from the soil and 16 reference sequences, while seven bacterial laccase fragments from the straw-amended soil and a reference sequence from giant panda feces were grouped into cluster D. Cluster G contained 12 bacterial laccase fragments associated with various genera of Proteobacteria, such as Pantoea, Klebsiella and so on, indicating high diversity of laccase-producing Proteobacteria in the straw-amended soil.

Cloning and expression of bacterial laccase genes from the straw-amended soil

The complete ORFs of five representative bacterial laccase genes were obtained using mTAIL-PCR (Table S2). Together with two full-length genes directly cloned from the soil metagenome, these seven bacterial laccases were translated into amino acid sequences. These seven bacterial laccases consisted of 445–629 amino acids, with molecular weights ranging from 49 kDa to 70 kDa and theoretical isoelectric points between 5.5 and 9.3. Among them, S_14142, S_63994, S_326900, and S_507105 contained TAT signal peptides, indicating their potential secretion via the twin-arginine translocation system (Table 3 and Fig. 5). Sequence alignments revealed low sequence identity among the seven laccases, with identities ranging from 23% to 52%, except for S_230270 and S_394411, which shared 70% identity (Table S3). Despite the differences at the whole protein level, the copper-binding sites were highly conserved (Fig. 5). None of the laccases exhibited more than 80% identity to previously reported sequences, confirming their novelty (Table S2).

For functional validation, S_63994, the most abundant gene, was designated lacS1 and expressed in E. coli without the signal peptide-coding sequence. Following microaerobic induction, recombinant lacS1 was purified using Ni²⁺-NTA metal-affinity chromatography, yielding a single protein band with an apparent molecular mass of approximately 48 kDa on 10% SDS-PAGE (Fig. S2). To verify the identity of the purified enzyme, the target band was excised and analyzed by LC-ESI-MS/MS, revealing eight peptides matching the recombinant lacS1 sequence and confirming its identity (Fig. S3).

For functional characterization of the recombinant lacS1, ABTS was used as the substrate. The enzyme was active across a broad pH range, spanning from 4.0 to 8.0, with optimal activity observed at pH 5.0 (Fig. 6A). It retained over 60% of its maximum activity after incubation at 37 °C for 1 h across pH values ranging from 5.0 to 9.0, demonstrating significant stability within this range (Fig. 6B). The optimal temperature for lacS1 activity was found to be 50 °C, and the enzyme exhibited thermostability at 40 °C (Figs. 6C–6D). However, lacS1 showed a gradual loss of activity when incubated at 50 °C for prolonged periods (Fig. 6D). The addition of copper ions to the reaction mixture had no significant effect on the enzyme’s activity (Fig. 7A; p < 0.05 by Tukey). In contrast, increasing concentrations of NaCl gradually reduced lacS1 activity, with the enzyme retaining over 42% activity at 250 mM NaCl. However, at concentrations exceeding 500 mM NaCl, the enzyme rapidly lost activity (Fig. 7B).