Integrative visual omics of the white-rot fungus Polyporus brumalis exposes the biotechnological potential of its oxidative enzymes for delignifying raw plant biomass

The genome of P. brumalis is enriched in genes coding for oxidative enzymes active on lignin

The draft genome sequence (94.9× coverage) of the monokaryotic strain BRFM 1820 is 45.72 Mb large and was assembled into 621 scaffolds and 1040 contigs, with a scaffold L50 of 0.36 Mb. In total, 18,244 protein coding genes were predicted. Expert annotation of the genes coding for enzymes active on lignocellulose revealed typical features of white-rot genomes, such as the presence of GH7 cellobiohydrolases (three genes) and AA9 LPMOs (17 genes) active on crystalline cellulose. The genome of P. brumalis also holds a commonly observed suite of enzymes active on hemicelluloses and is rich in enzymes active on acetylated xylooligosaccharides (CE16 acetyl esterases; 14 genes), and on pectin (two GH105 rhamnogalacturonyl hydrolase and 11 GH28 polygalacturonase coding genes). The number of oxidoreductases involved in lignocellulose degradation (AAs, 102 genes) is among the highest compared to sequenced genomes of the taxonomic order Polyporales [14]. We inspected whether CAZy gene family expansions occurred in P. brumalis using CAFE, a computational tool that provides statistical analysis of evolutionary changes in gene family size over a phylogenetic tree [15]. The results were consistent with a previous comparison of AA2 gene repertoires in Polyporales genomes, which showed a trend for larger AA2 gene families in the phlebioid clade that contains the deeply studied white-rot fungus Phanerochaete chrysosporium (16 AA2 coding genes; Fig. 1, [16]). However, we observed in P. brumalis an expansion of the AA2 gene family (19 genes, Additional file 1: Table S1), a feature also found in the core polyporoid clade for Trametes versicolor (26 genes) and Trametes cingulata (22 genes). AA2s are class II heme-containing peroxidases (PODs), i.e., MnPs, LiPs and VPs able to oxidize lignin [5] and non-ligninolytic Generic Peroxidases (GP: 1.11.1.7). Among them, 16 AA2s from P. brumalis were predicted to be secreted and could theoretically participate in the degradation of lignocellulosic substrates in the extracellular space. Remarkably, P. brumalis possesses a high number of VPs (nine genes) compared with other Polyporales species such as P. chrysosporium (which lacks VP genes), Pycnoporus cinnabarinus (two genes) and T. versicolor (two genes). VPs are characterized by wide substrate specificity due to the presence of three different catalytic sites in their molecular structure (i.e., a main heme access channel, a Mn²⁺ oxidation site, and a catalytic tryptophan exposed to the solvent) [17]. Their catalytic properties enable both the direct oxidation of phenolic and non-phenolic units of lignin, and Mn³⁺-mediated oxidations. Nine short MnPs complete the repertoire of ligninolytic AA2s. Unlike the long MnPs specific for Mn²⁺ (e.g., typical MnPs from P. chrysosporium), the members belonging to the short MnP subfamily can also oxidize phenols in Mn-independent reactions [18].

The AA2 gene family expansion in P. brumalis was accompanied by the expansion of the gene family AA3 (36 genes, excluding a unique cellobiose dehydrogenase coding gene that has the modular structure AA8-AA3). AA3 enzymes are flavo-oxidases with a flavin-adenine dinucleotide (FAD)-binding domain. Members of the family AA3 contribute to lignocellulose breakdown, among others through the generation of H₂O₂ that activates peroxidases [19] and the regeneration of AA9 LPMOs active on cellulose via the redox cycling of phenolics and the corresponding quinone mediators [20–22]. An additional role for AA3s could include contributing to the oxidative cleavage of polysaccharides through the generation of H₂O₂, which could subsequently be converted to hydroxyl radical OH^• by Fenton reaction as previously described in brown-rot fungi [19]. Twelve of these AA3s were predicted to be secreted and belong to the CAZy subfamily AA3_2 which includes predicted aryl alcohol oxidases (AAO, EC 1.1.3.7), aryl-alcohol quinone oxidoreductases (AAQO), glucose 1-oxidases (GOX, EC 1.1.3.4), glucose dehydrogenases (GDH, EC 1.1.5.9) and pyranose dehydrogenases (PDH; EC 1.1.99.29; [23–25]; Additional file 1: Table S1). There was no gene family expansion observed for AA9s in P. brumalis. To conclude, we observed the co-occurrence of gene family expansions for putatively secreted lignin-active peroxidases and auxiliary enzymes contributing to the generation of H₂O₂ and redox cycling of electron donors/acceptors, which could contribute to the oxidative degradation of wheat straw polymers by P. brumalis.

Polyporus brumalis selectively degrades lignin during solid state fermentation on wheat straw

The dikaryotic strain P. brumalis BRFM 985 was grown under solid-state cultivation with wheat straw as source of nutrient and support for 4, 10 or 15 days. Initial mycelium content in the packed bed bioreactor was of 5.8 mg (dw)/g substrate (dw). P. brumalis developed a matrix of hairy mycelia covering the wheat straw at Day 15 (Fig. 2a). The qPCR quantification of the fungal biomass (Fig. 2b) showed significant growth from Day 4 to Day 15 and an important weight loss (about 20%) of the total biomass within the bioreactors (Table 1). This result was consistent with previous work performed by Shi and co-workers [26]. As expected, wheat straw was selectively delignified by P. brumalis (Table 1). After 15 days culture, the ratios of lignin loss to cellulose loss (L/C) and lignin loss to holocellulose loss (L/H) reached 2.08 and 2.29, respectively. We investigated whether the changes in biochemical composition of wheat straw during fermentation had an impact on the accessibility of carbohydrates to glycoside hydrolases. This was evaluated using a commercial cocktail of enzymes. As shown in Table 1, the carbohydrate net conversion yields of fermented wheat straw increased over time of culture and reached 35% for cellulose at Day 15. For holocellulose, the carbohydrate net conversion yields were smaller (up to 30% at Day 15), which can be explained by the poor content of the cocktail in enzymes active on hemicelluloses.

	Weight loss (%)	Composition loss (%)			Selectivity		Carbohydrate net conversion yield (%)
		Cellulose	Holocellulose	Lignin	L/C	L/H	Cellulose	Holocellulose
Control	ND	ND	ND	ND	ND	ND	11.1 ± 0.1	10.5 ± 0.1
Day 4	0 ± 1.9	0 ± 1.8	0 ± 2.0	2.4 ± 1.8	ND	ND	19.6 ± 1.0	16.4 ± 0.6
Day 10	10.7 ± 0.1	10.1 ± 0.1	10.3 ± 0.1	18.6 ± 2.2	1.84	1.81	23.6 ± 0.4	20.9 ± 0.3
Day 15	19.6 ± 1.2	17.4 ± 1.2	15.8 ± 1.3	36.2 ± 1.0	2.08	2.29	35.2 ± 1.5	30.0 ± 0.3

ND not determined, L/C lignin to cellulose losses ratio, L/H lignin to holocellulose losses ratio

Polyporus brumalis activates an extensive arsenal of oxidative enzymes during growth on wheat straw

To identify enzymes responsible for the observed selective degradation of lignin, we investigated the evolution of transcriptomes and secretomes of P. brumalis over time by comparing a set of SSF on wheat straw for 4, 10, and 15 days. The robustness of the experimental setup was confirmed by the high consistency of the mean and distributions of the normalized log2 read counts of the transcripts among the biological replicates (Additional file 2: Figure S1). We constructed network maps for the fungal transcriptome and secretome using the visual multi-omics pipeline SHIN + GO, which allowed us to overview the genome-wide transcriptomic and secretomic activities and pinpoint molecular events of interest [7]. We made omics topographies (‘Tatami maps’) to: (1) visualize nodes containing genes with similar transcription patterns with the corresponding count of secreted proteins; (2) calculate the node-wise mean of the normalized transcript read counts in each condition; (3) identify gene clusters showing high transcription levels and high counts of secreted proteins at specific time points.

To analyze early adaptive responses of the fungus, we compared Tatami maps made from the SSF and liquid conditions. Triplicates of 10-day stationary liquid cultures on malt extract were compared with those from 4-day SSF on wheat straw with an initial 6-day liquid culture on malt extract. The change of the conditions from the liquid culture to SSF on wheat straw triggered rapid shifts in the transcriptome (Fig. 3). A total of 727 genes were differentially highly transcribed after 4-day SSF (Additional file 3: Table S2; Additional file 4: Dataset S2). According to KOG predictions [27], these genes were mainly related to metabolism, cellular processes and signaling, information storage and processing. The importance of intracellular signaling in this adaptive response was highlighted by the induction of 14 genes coding for predicted protein kinases and five for putative transcription factors. In addition, genes involved in detoxification and excretion were upregulated at the early time point including twenty, five, and four genes coding for predicted Cytochrome P450, Glutathione S-Transferases (GSTs), and ABC transporters, respectively, suggesting the active detoxification of compounds released from the biomass or newly generated during the wheat straw degradation (Additional file 5: Figure S2). The upregulation of 53 CAZyme coding genes at Day 4 was indicative of the early adaptive response of the fungus to solid-state cultivation and to wheat straw as a carbon source. There were three short MnPs and six VPs identified in the secretome at this early time point. The secretion of delignifying enzymes was accompanied with transcription up-regulation and secretion of enzymes active on cellulose, such as Lytic Polysaccharide Monooxygenases (five AA9 LPMOs identified in the secretome), endo- and exo-β-1,4 glucanases (one GH5_5, one GH5_22, one GH6 and three GH7s) and a GH131 β-1,3;1,4-endoglucanase. The coding genes of a CE1 (acetyl-esterase), a GH10 (xylanase), and a GH28 (polygalacturonase) were upregulated and the enzymes were present in the secretome, suggesting the concurrent depolymerization of hemicelluloses and pectin (Additional file 6: Table S3), despite pectin being present in very small quantities in wheat straws [28].

The time-course dynamics of the transcriptome and secretome during SSF were further investigated at Days 10 and 15. We made transcriptome-based Tatami maps at each time point showing the averaged transcription of the replicates (Fig. 4a) and the corresponding secretome-based Tatami maps displaying the count of secreted protein IDs detected per node (Fig. 4b). We observed a few groups of genes (5% of the analyzed protein coding genes) highly transcribed at Day 4 and later on downregulated, and global consistency between transcriptome maps at Days 10 and 15, reflecting that regulations related to the early adaptive response to cultivation under SSF were followed by a global stabilization of the transcription regulations for most of the nodes (Fig. 5). We also observed an apparent time lag in the intracellular processes between transcription and protein secretion. The correlation of the transcriptomes and secretomes showed that the transcriptomes at Day 4 were more positively correlated with the secretomes at Day 10 and 15. The same trend was seen for Day 10 transcriptome with Day 15 secretome (Table 2).

	Day 4 transcriptome	Day 10 transcriptome	Day 15 transcriptome
Day 4 secretome	0.53	0.48	0.44
Day 10 secretome	0.64	0.66	0.62
Day 15 secretome	0.64	0.67	0.64

Italics: Higher correlation coefficients observed between Day 4 vs Day 10/15 and Day 10 vs Day 15

Fungal strain and growth conditions

Polyporus brumalis BRFM 985 was obtained from the International Center of Microbial Resources (CIRM; https://www6.inra.fr/cirm). The fungus was maintained on 2% (w/v) malt extract, 2% (w/v) agar at 4 °C. Inoculum preparation and SSF procedures were performed as described by Zhou et al. [12]. Briefly, the fungal inoculum was prepared by grinding mats (Ultra-Turrax, 60 s, 9500 rpm) obtained from 7-day-old Roux flasks (malt extract 2%, w/v). Solid-state fermentations were performed in packed bed bioreactors (glass columns) on 20 g dw (dry weight) sterilized wheat straw (≈ 4 mm, Haussmann soft wheat, VIVESCIA, Reims, France) impregnated with 6 mg dw mycelium/g dw substrate, and 1 mL metal solution (CuSO₄, FeSO₄, and MnSO₄ 0.9 μmol each/g dw substrate). The glass-column reactors were incubated at 28 °C for 4, 10 or 15 days. Air stream was filtered, wetted and flow rate set to 0.5/v.v/min. Control wheat straw was incubated under the same conditions without inoculum. Assays were performed in triplicate.

Wheat straw biochemical characterization

Fermented and control wheat straw from three replicates were pooled and homogenized. Dry weights were measured from 1-g (wet weight) aliquots after drying overnight at 105 °C. Five grams (dw) aliquots were washed with deionized water (5% w dw/v) for 24 h at 4 °C, with shaking at 300 rpm, and filtered through GF/F filters (Whatman). The recovered solid extracts were freeze dried and stored at 4 °C for biochemical composition and enzymatic hydrolysis analyzes. The recovered water extracts were stored at − 20 °C for secretome analyzes. Biochemical composition (klason lignin, cellulose and hemicelluloses) of the fermented wheat straw was determined in duplicate according to the NERL method [42]. Carbohydrate accessibility to enzymatic degradation was evaluated by measuring carbohydrate net conversion yields as described by Zhou et al. [12]. Briefly, control or fermented wheat straw was subjected to mild alkali treatment, followed by a 72 h enzymatic hydrolysis with 12 FPU/g substrate (dw) of cellulases from Trichoderma reesei (GC220, Genencor Danisco) and 60 U/g substrate (dw) of β-glucosidase from Aspergillus niger (Novozyme SP188, Sigma). The released glucose and reducing sugars were quantified, respectively, using the Glucose RTU kit (Biomérieux) and the dinitrosalicylic acid method [43].

Fungal biomass quantification

Fungal biomass was quantified using a PCR-based method as described by Zhou et al. [44]. Briefly, genomic DNA was extracted from fermented wheat straw samples using the NucleoSpin PlantII kit (Macherey–Nagel, France) and quantified by qPCR amplification of a 150-bp fragment in the 5.8S conserved sequence. Standard curves were established using known quantities of genomic DNA from P. brumalis BRFM 985. The PCR cycle was as follows: 30 s at 95 °C, and then 5 s at 95 °C, 5 s at 60 °C for 39 cycles, followed by a melt curve step (65–95 °C, with 0.5 °C steps). The presence of a single amplicon was checked on the melting curve. All reactions were performed in triplicate, and all qPCR runs included a negative control without template. Quantification cycles (Cq) were determined using the regression mode of the Bio-Rad CFX Manager™ software (v 3.0).

Genome sequencing

To facilitate genome sequencing and assembly, the monokaryotic strain BRFM 1820 was obtained by de-dikaryotization [45] from the dikaryotic strain BRFM 985. The method for protoplast isolation was adapted from Alves et al. [46]. All steps were performed at 30 °C. Fifteen 7-day-old cultures (9 cm in diameter) of BRFM 985 grown on 2% (w/v) Malt extract, 2% (w/v) agar were recovered and homogenized for 20 s with a Waring blender in 100 mL of YM medium (10 g/L glucose, 5 g/L bactopeptone, 3 g/L yeast extract and 3 g/L malt broth extract). The homogenate was grown for 24 h with agitation at 130 rpm. The culture was again homogenized 20 s, diluted twice in YM, and grown for an additional 24 h with agitation. After washing with 0.5 M MgSO₄·7H₂O, the mycelium was harvested by centrifugation, weighted and resuspended in lysis buffer (0.5 M MgSO₄, 7H₂O, 0.03 M maleic acid pH 5.8 with 15 mg Caylase C4 per mL and per 0.25 g wet weight mycelium). The mycelium was incubated with gentle shaking (80 rpm) for 3 h and filtered through Miracloth. Protoplasts were collected by centrifugation at low speed and gently resuspended in 1 M sorbitol (10⁷/mL). Dilutions (10⁴/mL) were incubated overnight without shaking in regeneration medium (0.5 M MgSO₄,7H₂O, 20 g/L glucose, 0.46 g/L KH₂PO₄, 1 g/L K₂HPO₄, 2 g/L bactopeptone and 2 g/L yeast extract) and spread on MA2 medium. After culturing for 3–7 days, slow growing colonies obtained from the regenerated protoplasts were examined by microscopy to select the ones without clamp-connection. Fragments of mycelium were stained with DAPI to verify they were monokaryotic. From the 20 protoplast-derived monokaryons isolated, the monokaryotic line BRFM 1820, phenotypically similar to its dikaryotic parent, was selected for genome sequencing.

For genome, two (270 bp fragment and long 4 kb mate-pair) libraries were sequenced. For 270 bp fragments, 100 ng of DNA was sheared to 300 bp using the Covaris LE220 and size selected using SPRI beads (Beckman Coulter). For the 4 kb library, 5 µg of DNA was sheared using the Covaris g-TUBE™ (Covaris) and gel size selected for 4 kb. The sheared DNA was treated with end repair and ligated with biotinylated adapters containing loxP. The adapter ligated DNA fragments were circularized via recombination by a Cre excision reaction (NEB) and randomly sheared using the Covaris LE220 (Covaris). Sheared DNA fragments were processed for ligation to Illumina compatible adapters (IDT, Inc) using the KAPA-Illumina library creation kit (KAPA biosystems). For transcriptome, stranded cDNA libraries were generated using the Illumina Truseq Stranded RNA LT kit. Sequencing was performed on the Illumina HiSeq 2500 sequencer using HiSeq TruSeq SBS sequencing kits, v4, following a 2 × 150 bp (or 2 × 100 bp for 4 kb library) indexed run recipe. Each fastq file was QC filtered for artifact/process contamination. Genomic reads were assembled with AllPathsLG version R49403 [47]. RNA-Seq reads were assembled using Rnnotator v. 3.3.2 [48].

Gene functional annotation

The genome was annotated using the Joint Genome Institute (JGI) Annotation Pipeline and made publicly available via JGI fungal genome portal MycoCosm [49]. Proteins predicted to be secreted by SignalP 4.1. (threshold 0.34; [50]), TargetP 1.1 [51] or Phobius 1.01 [52] were analyzed with pscan to withdraw proteins targeted at the endoplasmic reticulum. The retrieved proteins with no transmembrane domain according to TMHMM Server 2.0. or a single transmembrane domain corresponding to the signal peptide were identified as predicted secreted. For CAZyme gene expert annotation, all putative proteins were compared to the entries in the CAZy database [2] using BLASTP. The proteins with E values smaller than 0.1 were further screened by a combination of BLAST searches against individual protein modules belonging to the AA, GH, GT (Glycosyl Transferases), PL, CE and CBM classes (https://www.CAZy.org/). HMMer 3 [53] was used to query against a collection of custom-made hidden Markov model (HMM) profiles constructed for each CAZy family. All identified proteins were manually curated. Within families, subfamilies were manually defined according to their homology relationships between members of the focal family. Class II ligninolytic peroxidases (AA2s) were annotated as LiP, MnP or VP on the basis of the presence or absence of specific amino acid residues at the substrate oxidation sites [17] after homology modeling using crystal structures of related peroxidases as templates and programs implemented by the automated protein homology modeling server “SWISS-MODEL” [54].

Proteins classified as KOG “Cytochrome P450” were subject to BLASTP comparisons against the Fungal CYP Database (http://p450.riceblast.snu.ac.kr, last accessed June 10, 2018) [55]. These predicted CYPs were assigned to the corresponding family types based on their highest sequence similarity (at least 40%) against all fungal CYPs as followed by the International P450 Nomenclature Committee.

Gene family expansions/contractions

CAZyme gene family expansions or contractions were analyzed using the CAFE Software v3.1 (Computational Analysis of gene Family Evolution; [15]) and the numbers of genes coding for CAZymes from 18 Polyporales genomes. The selected public genomes were those from Pycnoporus cinnabarinus [56], Phanerochaete chrysosporium [57], Fomitopsis pinicola, Wolfiporia cocos and Trametes versicolor [58]. Genomes from Artolenzites elegans BRFM 1663 and BRFM 1122, Irpex lacteus CCBAS Fr. 238 617/93, Leiotrametes sp. BRFM 1775, Pycnoporus coccineus BRFM 310 and BRFM 1662, P. puniceus BRFM 1868, P. sanguineus BRFM 1264, Trametopsis cervina BRFM 1824, Trametes cingulata BRFM 1805, T. gibbosa BRFM 1770 and T. ljubarskii BRFM 1659 were newly sequenced at the JGI. The Russulales Stereum hirsutum [58] and Heterobasidion annosum [59] were used as outgroups. The protein sequences deduced from these genomes were downloaded from the Mycocosm portal (https://genome.jgi.doe.gov/programs/fungi/index.jsf), with authorization of the principal investigators when not yet published. Groups of orthologs were obtained from these genomes using orthoMCL version 2.0.9 [60] with default parameters and using NCBI BLAST version 2.4.0 + [61] in conjunction with MCL version 14-137 (http://micans.org/mcl/) under default settings. A set of 20 uni-copy genes was selected among the 25 best performing genes for Polyporales phylogenetics ([62]; Additional file 9: Dataset S1). The concatenated 20-gene dataset in each genome had a mean number of 28,996,9 amino acids, the shortest being for Pycci1 (28144 aa) and the longest for Stehi1 (30313 aa). Alignments were done using Mafft v7.271 [63] and filtered with GBLOCKS version 0.91b [64]. The final tree was generated with RaxML version 8.2.4 [65] and PROTGAMMAWAG as substitution model.

Identification of the fungal proteins secreted during fermentation

The water extracts collected after 4, 10 and 15 days of culture were filtered using 0.22-μm pore-size polyethersulfone membranes (Vivaspin, Sartorius), diafiltered with 50 mM sodium acetate (pH 5.0) and concentrated using a Vivaspin polyethersulfone membrane with a 10-kDa cutoff (Sartorius). LC–MS/MS analysis of the secretomes was performed as described by Navarro et al. [66]. Briefly, 10 μg of proteins was in-gel digested using a standard trypsinolysis protocol. For protein identification, online analysis of peptides was performed with a Q-exactive mass spectrometer (Thermo Fisher Scientific), using a nanoelectrospray ion source. MS/MS data were queried against the catalog of predicted proteins from the P. brumalis genome, and an in-house contaminant database, using the X!Tandem software (X!Tandem Cyclone, Jouy en Josas, France). Peptides that matched with an E value < 0.05 were parsed with the X!Tandem pipeline software. Proteins identified with at least two unique peptides and a log (E value) < − 2.6 were validated.

Transcriptome and secretome network maps

Total RNA was extracted from 100 mg tissue as described in Couturier et al. [67]. RNA quantity and quality were determined using the Experion RNA Std-Sens kit (QIAGEN). RNASeq was done using 150 bp long paired end reads obtained on Illumina HiSeq-2500 (Beckman Coulter Genomics). The quality of raw fastq reads was checked with FastQC (Simon Andrews, Babraham Bioinformatics, 2011; https://www.bioinformatics.babraham.ac.uk/projects/) and reads were cleaned with Sickle Tool (Joshi and Fass, 2011; https://github.com/najoshi/sickle) with the following criteria: quality threshold of 20; length threshold of 20. RNA reads were aligned to the genome of P. brumalis BRFM 1820 using TopHat2 with only unique mapping allowed. Read counts were determined by HTSeq, normalized using ddsNorm from the DESeq 2 Bioconductor package [68] and log2 transformed.

Genome-wide profiles with integrated transcriptome and secretome data were constructed using the SHIN part of the pipeline SHIN + GO [7]. A single self-organizing map (SOM) was generated with the normalized log2 read count of the genes from all biological replicates. The map units used were 437. The number of iterations used was 43,700 (437 × 100) times to minimize the Euclidian distances between the nodes for the optimal convergence. Two sets of omics topographies (Tatami maps) were generated based on the master SOM. Transcriptomic maps showing node-wise mean transcription were created with the averaged read count of the genes clustered in each node. Secretomic maps were made with the count of detected proteins. For each protein ID, the presence of the protein in the secretome was counted as one. The total count of secreted proteins was calculated node-wise according to the trained SOM. A summary map indicating condition-specific highly transcribed genes was also made. Genes in the nodes with a mean read count > 12 log2 were considered as highly transcribed in each condition at each time point. Such genes constituted approximately above 75th percentile of the entire gene population.

Spearman’s rank correlations of fungal transcriptome and secretome were estimated at individual time points using the node-wise mean transcription and the node-wise count of secreted proteins identified. Correlations were calculated for: (1) the entire genes; and (2) the subset of genes coding for predicted secreted proteins.