Wood modification by furfuryl alcohol caused delayed decomposition response in Rhodonia (Postia) placenta

The aim of this study was to investigate differential expression profiles of the brown rot fungus Rhodonia placenta (previously Postia placenta) harvested at several time points when grown on Pinus radiata (radiata pine) and P. radiata with three different levels of modification by furfuryl alcohol, an environmentally benign commercial wood protection system. For the first time the entire gene expression pattern of a decay fungus is followed in untreated and modified wood from initial to advanced stages of decay. Results support the current model of a two-step decay mechanism, with an initial oxidative depolymerization followed by hydrolysis of cell-wall polysaccharides. The wood decay process is finished, and the fungus goes into starvation mode after five weeks when grown on unmodified P. radiata wood. The pattern of repression of oxidative processes and oxalate synthesis found in P. radiata at later stages of decay is not mirrored for the high furfurylation treatment. The high treatment level provided a more unpredictable expression pattern throughout the entire incubation period. Furfurylation does not seem to directly influence the expression of core plant cell wall hydrolyzing enzymes, as a delayed and prolonged, but similar pattern was observed in the P. radiata and the modified experiments. This indicates that the fungus starts a common decay process in the modified wood, but proceeds at a slower pace as access to the plant cell wall polysaccharides is restricted. This is further supported by the downregulation of hydrolytic enzymes for the high treatment level at the last harvest point (mass loss 14%). Moreover, the mass loss does not increase the last weeks. Collectively, this indicates a potential threshold for lower mass loss for highly modified wood. IMPORTANCE Fungi are important decomposers of woody biomass in natural habitats. Investigation of the mechanisms employed by decay fungi in their attempt to degrade wood is important for both the basic scientific understanding of ecology and carbon cycling in nature, and for applied uses of woody materials. For wooden building materials long service life and carbon storage is essential, but decay fungi are responsible for massive losses of wood in service. Thus, optimizing durable wood products for the future are of major importance. In this study we have investigated the fungal genetic response to furfurylated wood, a commercial environmentally benign wood modification approach, that improves service life of wood in outdoor applications. Our results show that there is a delayed wood decay by the fungus as a response to furfurylated wood and new knowledge about the mechanisms behind the delay is provided.


138
In the following experiments the P. radiata wood was modified at three different levels. 139 The three modification levels had a mean Weight Percent Gain (WPG) of 3.8±0.7%, 24.0±3.5%, 140 36.6±5.0%. For simplicity the treatments are named WPG4, WPG24 and WPG37, respectively.

141
The experiment on the unmodified wood is named P. radiata. 142 Mass loss calculations. Five weeks after inoculation of R. placenta strain FPRL 280 on the 143 P. radiata wood blocks, they had a mean mass + treatment loss of 28.8±4.0%, while the modified 144 WPG37 only had a mean mass + treatment loss of 13.5±3.3% after 21 weeks (Fig. 1). Mass loss is 145 within the wood protection literature mostly referred to as the mass loss of the entire wood + 146 treatment system, i.e. dry weight of the treated wood before fungal inoculation compared to dry 147 weight after fungal decay. Another way to measure mass loss is to assume that only the wood is 148 decayed and eliminate the WPG of the treatment from the mass loss calculation. Both methods 149 resulted in the same general trends, but the exclusion of WPG results in slightly higher mass losses 150 as expected (Fig. S1). The differences between the two mass loss calculation approaches for the 151 three furfurylation levels at the last harvesting points were (wood + treatments vs. only wood): 152 WPG4 week 9 41.7% vs. 42.5%, WPG24 week 21 29.0% vs. 34.6%, WPG37 week 21 13.5% vs. 153 16.0%. We use the wood + treatment results as presented in Fig. 1 as our mass loss detection 154 system.
155 Differential gene expression. In order to investigate the genetic basis of the different 156 behavior of R. placenta growing on different levels of furfurylated wood we sequenced 1, 2, 3, 4 and 5, WPG4 -weeks 3, 6 and 9, 6,9,12,15,and 18, weeks 3, 6, 9, 12, 15 and 18. 161 Even if the genome and transcriptome of the American R. placenta strain MAD 698-R has 162 been sequenced, this could not be used for mapping purposes in this study. The genome sequenced 163 American strain MAD 698-R and the European strain FPRL 280 used in this study are significantly 164 different, with a mapping success of only 40-60% when the reads were mapped to the genome of 165 MAD 698-R (results not shown). We therefore produced Illumina Hiseq paired-end data to be used 166 for a transcriptome assembly. The resulting assembly had 56 520 contigs (transcripts), and covered 167 99.3% of the conserved BUSCO fungal genes (Table S1; Table S2). In the transcriptome, 9 355 168 contigs had an annotation from Blastx, and 18 917 contigs were given a PFAM annotation. In all 169 further analyses the sequence data were mapped to this annotated transcriptome assembly. The 170 mapping success of sequence data to this transcriptome was more than 90% for all libraries.

171
The gene expression data demonstrate consistent results for each treatment with a gradient 172 clustering of replicate samples according to treatment with few outliers (Fig. 2). There is some 173 biological variation between replicates as expected from natural variation in wood and the variation 174 in modification of the wood blocks. The P. radiata and WPG37 experiments display a wider range 175 of biological variation compared WPG24 and WPG4. i.e. WPG4 vs P. radiata, WPG24 vs P. radiata and WPG37 vs P. radiata. (Table 1). Functional enrichment analyses of the resulting gene lists were inferred using the annotated PFAM domains 182 and GO terms from the annotated transcriptome. Compared to P. radiata, the WPG4 showed 183 upregulation of zinc-binding dehydrogenase domains (PF00107.21) and two GO terms related to 184 zinc ion binding and oxidoreductase activity. The same enrichment was also found upregulated in 185 WPG24 and WPG37 when compared to P. radiata. In addition, these two latter treatments showed 186 upregulation of more oxidoreductase domains, most in WPG37. No PFAM domains or GO terms 187 were found downregulated in WPG4 or WPG24. However, WPG37 showed a strong 188 downregulation of eight PFAM domains and seven GO terms with functions related to protein and 189 peptide degradation, the ubiquitin-proteasome pathway and the Ras gene family compared to P. 190 radiata (Table 1), e.g. two domains related to proteasome (PF10584.4 and PF00227.21) and Ras 191 (PF08477.8 and PF00071.17). These terms were not found among the other comparisons.

192
As an alternative analysis method to the multifactorial DE analyses above, we also clustered 193 the genes with similar expression profiles. The read counts were grouped into ten clusters (K) for 194 each treatment, and samples from P. radiata, WPG4, WPG24 and WPG37 were analyzed 195 separately (Fig. 3). For these clusters we did functional enrichment with PFAM and GO terms 196 (Table S3), and also investigated the placement of known genes related to plant cell wall decay in 197 these clusters (Table S4).

198
For most of the clusters, an observed pattern of higher expression in one or a few of the 199 harvest points was found (Fig. 3). For P. radiata, two clusters were directly linked to wood decay 200 and carbohydrate active enzymes (Table S3). One of these clusters (P. radiata-K2) was related to 201 early depolymerization of hemicellulose and pectin (enriched for GH28 and GH43 domains and 202 containing the genes OxaD, Man5a and CE16b; see Table 2 for information about abbreviated gene 203 names) and highly expressed in week 1, while the other cluster (P. radiata-K5) was related to later 204 stages of cellulose depolymerization (enriched for GH3 and hydrolase activity, and containing the 205 genes Cel5b, Cel2, bGlu, Xyl10a, bXyl; Table 2) and was highly expressed in week 2. Similar 206 enrichments were found for two clusters of WPG4 (WPG4-K2 which is highly expressed early, 207 containing Man5a, CE16a and Gal28a and WPG4-K5 which was expressed late in contrast to P. 208 radiata and containing OxaD and Xyl10a), however no clusters can be directly compared across 209 treatments. Two clusters in WPG24 also have enrichment for GO terms related to hydrolase 210 activity; WPG24-K3, which is highly induced in week 3 and contains Gal28a and CE16b and 211 WPG24-K9, which is a large group with no specific induction pattern across harvesting points, 212 where all the other specific genes we investigated were placed. The same pattern was demonstrated 213 for WPG37, where all specific genes, except CE16b, were placed in a large group (WPG37-K10) 214 with no specific induction time. WPG37-10 also had enrichment for GO terms such as protein 215 binding, oxidation-reduction processes and catalytic activity (see Table S3 for more details).

216
Functional enrichment was only found for one other cluster of WPG37 (WPG37-K1). This cluster 217 was enriched for functions related to salt and water stress and was highly expressed in the second 218 harvest point, week 6. One cluster with one or a few of the same domains as in cluster WPG37-K1 219 could be found for all the other treatments (P. radiata-K1, WPG4-K6 and WPG24-K3). For the 220 WPG24 this was induced in the first harvest point, thus earlier than for WPG37. In P. radiata-K1 221 and WPG4-K6 there were no clear pattern of induction time for clusters with these stress domains 222 ( Fig. 3; Table S3).

223
Differential gene expression and functional enrichment within treatments. Pair-wise 224 comparisons of different harvest points within treatments revealed large gene expression 225 differences between the first and the later harvest points, for all treatments. In the P. radiata experiment, 2450 R. placenta genes were differential expressed between week 1 and week 2, while 227 only 10 genes were differential expressed between week 2 and week 3 (Table S5).

228
The furfurylated wood modification demonstrated a similar temporal pattern of DE genes 229 as the P. radiata experiment (Table S6). For WPG4 there were 103 DE genes between week 3 and 230 week 6 and 481 DE genes between week 3 and week 9. No DE genes were observed between week 231 6 and week 9. For WPG24, large numbers of DE genes were observed between week 3 and the 232 later time points, but again very few in between the later harvest points. For WPG37, few genes 233 were differentially expressed (between harvest points) in general. However also here, more DE 234 genes were found between week 3 and the later harvest points, and fewer among the later harvest 235 points.

236
Functional enrichment analyses of the pair-wise analyses conformed with the expected two-237 step decay in this system (Table S7). These results demonstrated that for P. radiata experiment 238 there was an increased expression of functions related to hydrolase and catalytic activity from week 239 1 to the later harvest points, especially pronounced in week 2 and week 3 (Table S7). This was also 240 found, but to a lesser extent in WPG4 and WPG24. In WPG4, most of this response was found 241 between week 3 and week 9, while for WPG24, the response mainly started in the comparison 242 between week 3 and week 12 to week 18. For WPG37, this was not pronounced, and no hydrolase 243 activity was enriched. In the first harvest points of the P. radiata, WPG4 and WPG24 treatments 244 we found less enrichment of upregulated functions with the exception of some upregulation of GO 245 terms related to protein binding and metal ion binding. However, in WPG37 there was an 246 enrichment of several terms related to iron binding, heme binding, oxidoreductase activity and 247 cytochrome P450s in the upregulated gene set of week 3 compared to week 18. Significant enrichment of a GO term of oxidation-reduction process was found in later stages in WPG4 and 249 WPG24.

250
As week 1 in P. radiata was considered the first initial step of decay, with only 0.8% mass 251 loss, we compared P. radiata week 1 to all other time points and treatments (Table S8). As with 252 the timeseries analyses, a signal of downregulation of protein degradation was found in all the 253 WPG37 harvest points compared to P. radiata. This was seen as enrichment of protein kinases, 254 proteasome, F-box domain and endopeptidase activity in P. radiata. This pattern was also found 255 in the early harvest points of WPG24. In contrast week 1 of the P. radiata showed upregulation of 256 wood decay related functions as sugar transporters, GHs, and dehydrogenases.  Detailed qRT-PCR analyses of key genes of interest. Large RNAseq datasets tend to 267 contain a lot of biological variation as seen in our PCA plot, which will affect downstream 268 differential expression analyses. qRT-PCR serves both as a RNAseq control but also provides 269 deeper/more detailed insight into the expression of individual genes. Here we have used qRT-PCR primers for key genes involved in wood decay and applied a traditional qRT-PCR approach (Thus 271 the same genes as those reported in the cluster analyses; Table 2; Table S4; Table S9). Plots of both 272 RNAseq and qRT-PCR data for selected genes are provided ( Fig. S3-S7). The results below are 273 based on these qRT-PCR data, RNAseq data are commented on when the trends deviated from the   Table   300 2), two copper radical oxidases (Cro1, Cro2; Table 2) and a benzoquinone reductase (BqR; Table   301 2).

302
GMC oxidoreductases are flavoenzymes that oxidize a wide variety of alcohols and 303 carbohydrates, with concomitant production of hydrogen peroxide or hydroquinones (36). AOx1 304 and AOx2 were in the present study upregulated in WPG37 compared to P. radiata, and for AOx2 305 WPG24 was also upregulated. AOx3 was highly expressed but no significant differences in 306 expression were found within (except an upregulation WPG24 week 3) or between treatments.

308
Copper radical oxidases (CAZy family AA5) are widely distributed among wood decay 309 fungi (37), and oxidize a variety of substrates and produce H2O2. Cro1 and Cro2 were in the present 310 study upregulated week 1 for P. radiata and week 21 for WPG24. Since RNAseq was not 311 performed beyond week 18 this trend was not confirmed. Cro2 was upregulated in WPG24 312 compared to WPG4.

313
Benzoquinone reductases (CAZy family AA6) are intracellular enzymes suggested to have 314 a role in degrading lignin and derived compounds in white rot fungi but also suggested to provide protection from oxidative stress by acting as redox active toxins and also contribute to oxidative 316 depolymerization of wood cell wall polymers via the production of hydroquinone chelators (38).

317
In the present study BqR was upregulated in later stages of decay both for P. radiata and for 318 WPG37.

319
Hydrolytic enzymes involved in polysaccharide depolymerization: Hemicellulose and 320 pectin degradation. Fig. S5 illustrates the selected genes involved in hemicellulose hydrolysis.

321
One endomannanase in CAZy family GH5 (Man5a; Table 2) was found to be upregulated in P. 322 radiata versus WPG37. All treatments showed downregulation with increasing incubation time.

323
For the two endoxylanases in CAZy family GH10 (Xyl10a and Xyl10b; Table 2) there was 324 a tendency for upregulation in furfurylated wood at 3 weeks incubation. For RNAseq an 325 upregulation in P. radiata week 2.

330
Carbohydrate esterases catalyze the de-O or de-N-acylation of substituted saccharides. The 331 selected carbohydrate esterase CE16a (Table 2) was upregulated in P. radiata compared to WPG4 332 and WPG37. All treatments showed downregulation with increasing incubation time. The 333 carbohydrate esterase CE16b (Table 2) provided very low expression levels with qRT-PCR 334 (especially P. radiata and WPG4). A clearer trend was found with RNAseq, especially for WPG24 335 and WPG37 with decreasing expression levels with increasing incubation time.

344
Glucoside hydrolase XyGEg (Table 2) shows high sequence similarity to several 345 endoglucanases with activity on cellulose and xyloglucan. This enzyme was upregulated in WPG4 346 vs. WPG37. In P. radiata samples it was a gradual increase in XyGEg with increasing incubation 347 time.

Expansins with a possible role in loosening plant cell-wall interactions.
Expansins are 357 hypothesized to increase enzyme accessibility by loosening plant cell-wall interactions, although 358 no catalytic mechanism is known (39-41). The two selected genes predicted to encode expansins 359 (Exp1 and Exp2; Table 2) are given in Fig. S7. The expression patterns of the two expansins did 360 not show a clear trend, but for Exp1 in WPG37 the expression was upregulated week 3 and 9 361 compared to later decay stages. 364 We have investigated if wood modification by furfuryl alcohol causes altered decomposition 365 response in Rhodonia (Postia) placenta, using RNAseq and qRT-PCR, in unmodified and 366 furfurylated P. radiata wood.

367
From our four treatments (P. radiata, and the three levels of modification with furfuryl 368 alcohol polymer) we found that WPG4 closely followed the trend on P. radiata, while WPG24 369 seemed to be an accelerated version of the decay pattern of WPG37. For simplicity, we mainly 370 focus on P. radiata decay processes and the comparison of P. radiata decay vs. WPG37 decay in 371 the discussion.

372
In our P. radiata decay experiments we observed wood decay from incipient stages with 373 mean mass loss of 0.8% in week 1 to severely decayed wood with a mean mass loss of 29% at 374 week 5. Our results support the previously suggested two-step decay mechanism of brown rot fungi 375 (13,16,17,19). Week 1 of P. radiata demonstrated an early decay transcriptome response, with 376 higher expression of oxidative enzymes, polygalacturanase GH28 and oxalate synthesis genes in 377 both the RNAseq data and for the specific genes subject to qRT-PCR ( Fig. 4; Fig. S4).

378
This early response was followed by a polysaccharide depolymerization stage, with 379 upregulation of hemicellulases, e.g. endoxylanases (Xyl10b and Xyl10a), beta xylosidase (bXyl) 380 and endomannanase (Man5a) in week 2 and week 3. In our study there was an abrupt 381 downregulation of many core hemicellulases in week 4 during the later stages of decay (Fig. 4, Fig.   382 S6); Xyl10a, Xyl10b, bXyl and Man5a). This is in agreement with observations made by Zhang 383 and Schilling (19). It has been demonstrated that soluble sugars, in particular cellobiose, acts as the 384 main inducing agent in the switch from oxidative (week 1) to hydrolytic depolymerization (19, 42).
Cellulose active enzymes were also induced in week 2, and were further upregulated in 386 week 3. This was particularly obvious from the qRT-PCR expression pattern of the GH3 387 betaglucosidase (bGlu), and the LPMO (Fig. S6). Some cellulases are also highly expressed in the 388 more advanced stages of decay, e.g. both the endoglucanase GH5 (Cel5b) and the glucoside 389 hydrolase GH12 (Cel12a) are highest expressed in week 5. This high expression of cellulases at 390 late decay stage was also observed by Zhang and Schilling (19). Thus, when the more easily 391 accessible hemicelluloses have been degraded, expression of cellulases Cel5b and Cel12a are 392 maintained at a high level, well into the phase where the fungus is experiencing starvation 393 (discussed further down).

394
When the decay of the untreated P. radiata was compared to the modified wood, we 395 observed that the oxidative enzymes were highly expressed in the first harvest point (week 3) also 396 for the modified WPG37. This is confirming the finding in Alfredsen et al. (28) where R. placenta 397 expressed high levels of oxidative enzymes and produced oxalate during 8 weeks colonization of 398 furfurylated Scots pine. However, the pattern of repression of these oxidative processes and the 399 oxalate synthesis found in P. radiata at later stages is less clear (Fig. 4). This is particularly obvious 400 when comparing the expression of oxalate decarboxylase (OxaD) and AA3 GMC oxidoreductase 401 (AOx1), which showed a distinct time of induction on P. radiata but not on WPG37 where elevated 402 expression was randomly distributed over time between samples. This can conceivably be 403 explained by a repeated exposure to regions with heavily modified substrate, thus re-inducing 404 oxidative processes to overcome furfurylation in order to expose degradable substrate to the fungus 405 allowing further growth, or a direct response to the furfuryl polymer.  Thus, based on the observation that a delayed but, similar pattern was observed in both the P.

439
As a support for the notion that the furfurylation is non-toxic for the fungus, we found no 440 evidence for expression of more defense mechanisms in the modified wood. Rhodonia placenta is 441 known to be highly tolerant to substances as copper, mainly due to the ability to produce oxalic 442 acid that chelate and precipitate copper and other transition metals (49). However, in addition to 443 these more general functions that are difficult to separate from wood decay mechanisms, copper 444 tolerant fungi are also known to express catalases and ATP-pumps related to copper transport in 445 response to toxins (50). These functional categories were not enriched in the differential expressed 446 gene sets in our study when comparing WPG37 to the P. radiata, indicating that the fungus is not 447 experiencing a more toxic atmosphere in the modified wood.

448
Notably, the most pronounced transcriptome differences in the P. radiata compared to the 449 modified wood is the strong induction of functions related to the ubiquitin/proteasome pathway, 450 protein degradation and the RAS pathway. All these functions are related to carbon starvation and 451 are highly enriched in the late decay of the unmodified P. radiata. Thus, we hypothesize that R. 452 placenta growing on P. radiata is starving and that the fungus has consumed the majority of available carbon sources in the wood in the latest harvest points of the P. radiata experiment. This 454 is further supported by the downregulation of CAZymes the last two weeks in this experiment.

455
The ubiquitin/proteasome pathway is a conserved pathway in all eukaryotic organisms, and 456 is important to various cellular processes as recycling of intracellular protein (where unnecessary 457 proteins are degraded to amino acids that can be reused to produce new proteins) and programmed 458 cell death. Recently, the pathway was shown to be activated by carbon starvation in the 459 ectomycorrhizal basidiomycete species Paxillus involutus (51). In P. involutus, 45% of the 460 transcripts were differentially regulated during carbon starvation. This large response is also shown 461 in our study where most enriched functions in P. radiata compared to the WPG37 for the time 462 series analyses could be connected to the U/P or the Ras pathway. This could reflect the higher 463 need for translocation of nitrogen rich resources in the P. radiata substrate to areas that need new 464 protein synthesis than is the case in WPG37 or simply reflect a general starving response due to 465 lack of substrate.

466
The other pattern observed as upregulated in P. radiata compared to WPG37 was several 467 domains related to Ras proteins. These proteins are involved in cell proliferation and growth of the 468 mycelia. Carbon starvation supports an upregulation of these proteins. In starving mycelia it has 469 been shown that the diameter of the hypha is reduced, while it grows to cover larger area (51, 52). front. This wafer approach works relatively well to separate different initial decay stages. For more 488 advanced stages of decay a different approach is needed. In the current study we selected small and 489 homogeneous (earlywood) samples in order to enable fast colonization and as homogenous solid 490 wood substrate for decay as possible. Still, some variation of gene expression is expected as new 491 colonization pathways will be found close to areas with previously colonized wood. This effect is 492 expected to increase with increasing WPG level, since areas long colonized will cease to provide 493 nutrients able sustain survival while the others more recently invaded will provide nutritious 494 substrates. In our experiments, the wood modification treatment modifies the composition of the 495 wood, and the fungus is therefore forced to respond differently than when it encounters unmodified wood. Hence, comparison at similar mass losses between treatments was not the goal of this study, 497 but rather the shift of gene expression over time between and within the different treatments.

498
Increased knowledge about brown rot decay mechanisms is important for an expanded 499 understanding of the fungal decay process in general since ecosystem carbon flows are closely 500 linked to wood degrading fungi. From an industry perspective the findings in this study show that 501 successful inhibition of the initial oxidative decay is a clue to success for future wood protection 502 systems.

504
This is the first time the entire gene expression pattern of a decay fungus is followed in untreated 505 and modified wood from initial to advanced stages of decay. From these observations, we have 506 demonstrated how the furfurlyated modification delays the fungus gene expression while growing 507 on this substrate. All treatments (modified and unmodified) expressed the expected staggered decay 508 mechanism, i.e. oxidative enzymes earlier, and then the hydrolyzing enzymes later. Thus, we show 509 that the fungus experiences a similar substrate situation on the modified wood, or at least is 510 following the same modality of gene expression. The major difference is that the responses are 511 delayed and elongated compared to unmodified wood. However, from the downregulation of 512 CAZymes in the latest harvest points in the WPG37, the fungus seems to be finalizing the decay 513 also in the furfurylated wood at a mass loss of only 14%, compared to 29% in P. radiata. This 514 suggests that the fungus never get access to the remaining carbohydrates from the modified wood 515 even at low mass loss. The lower levels of modification in the WPG4 and WPG24 treatments show 516 intermediate expression differences and higher mass loss. The mode of action of the modification 517 is still uncertain. However, the similar decay mechanisms, and the lack of expression of defense mechanisms may indicate that the furfurylation functions as both a physical barrier and a factor 519 that creates a less hydrated environment for the fungus. These hypotheses should be investigated 520 to further improve environmentally friendly modification processes.

523
Wood material. In order to get as homogeneous wood material as possible plugs (Ø = 6 524 mm, h = 10 mm ) were prepared from Pinus radiata D. Don earlywood according to Beck et al. and initial dry mass after treatment the samples were dried at 103 °C for 18 hours then cooled down 535 in a desiccator before the dry weight was recorded. The samples were left in a climate chamber at 536 65% relative humidity and 20 °C until stable weight before they were wrapped in sealed plastic 537 bags and sterilized by gamma irradiation (25 kGY directly into a container with liquid nitrogen. The samples were then stored at -80°C.

568
They were then cooled in liquid nitrogen again before a second round of grinding.  sequencing.uio.no). All samples were prepared with the strand specific TruSeq TM RNA-seq library 590 preparation method. In order to produce a high quality de novo transcriptome we sequenced the 591 strain grown on a 4% (w/v) Difco TM malt agar (VWR) media in addition to ten selected samples 592 spanning the entire experimental setup on one lane on the Hiseq2500 producing 125 bp paired end 593 sequences. These libraries generated ~395 million read pairs, which after quality control and 594 trimming was reduced to ~250 million read pairs.  (Table S1). 615 We also evaluated the assembly using BUSCO v. 2.0 using the recommended parameters  (Table S2).

618
Assembly annotation using Trinotate and the JGI Fungi Portal. The assembly was 619 subjected to the Trinotate annotation pipeline according to the manual which provided a generic 620 annotation of assembly. However, as there is a study (17) looking into key wood-degrading genes 621 using the annotation from the R. placenta genome, we manually annotated our assembly using the to explore that overall data and to plot raw counts.

638
The overall workflow with code is available in the supplementary information, but briefly 639 we followed the steps described below:

640
The coldata object describing the overall experiment was set up with time, treatment, mass 641 loss and condition where condition in practice is the intercept between time and treatment. We the data using the extra term(s) in the full model is more than expected.
In edgeR we opted for a ~0 + treatment+time formula without generating an intercept term.

653
When estimating overall dispersion, we used the robust=TRUE option to better handle outliers in 654 the data. For the multifactorial test itself we chose to replace the standard glmFit with glmQLFit 655 which uses a quasi-likelihood F-test on the likelihood ratio statistics instead of the chisquare 656 approximation. In this way we should obtain a more conservative control of the type I error rate as 657 it takes into account the uncertainty in estimating dispersion for each gene -especially when the 658 number of replicates is small.

659
To enable exploration of specific contrasts between given conditions we also ran a ~0 + 660 condition formula in edgeR using the same setup as described above.

661
We evaluated each treatment vs the P. radiata over time from the treatment+time 662 multifactorial analysis. For the pairwise contrasts between harvest time within treatments we used 663 the ~condition multifactorial analysis.

664
The overall exploration of the data revealed that the WPG24 treatment series has a few 665 outliers. RNAseq data is highly variable and the individual variation between replicates can be 666 large. We found that 2 replicates in the WPG24 3 weeks condition were relatively deviating from 667 the other two replicates. We reran the above described analysis without the most extreme outlier 668 (WPG24 3 weeks replicate 2) and observed a large change in the number of reported significant 669 differentially expressed genes. However, as these differences not necessary improved the results

670
we decided to keep all samples in the final dataset.