Weighted Gene Co-expression Network Analysis Identifies Critical Genes for the Production of Cellulase and Xylanase in Penicillium oxalicum

Genes involved in cellular processes undergo environment-dependent co-regulation, but the co-expression patterns of fungal cellulase and xylanase-encoding genes remain unclear. Here, we identified two novel carbon sources, methylcellulose and 2-hydroxyethyl cellulose, which efficiently induced the secretion of cellulases and xylanases in Penicillium oxalicum. Comparative transcriptomic analyses identified carbon source-specific transcriptional patterns, mainly including major cellulase and xylanase-encoding genes, genes involved in glycolysis/gluconeogenesis and the tricarboxylic acid cycle, and genes encoding transcription factors, transporters and G protein-coupled receptors. Moreover, the weighted correlation network analysis of time-course transcriptomes, generated 17 highly connected modules. Module MEivory, comprising 120 members, included major cellulase and xylanase-encoding genes, genes encoding the key regulators PoxClrB and PoxXlnR, and a cellodextrin transporter POX06051/CdtC, which were tightly correlated with the filter-paper cellulase, carboxymethylcellulase and xylanase activities in P. oxalicum. An expression kinetic analysis indicated that members in MEivory were activated integrally by carbon sources, but their expressional levels were carbon source- and/or induction duration-dependent. Three uncharacterized regulatory genes in MEivory were identified, which regulate the production of cellulases and xylanases in P. oxalicum. These findings provide insights into the mechanisms associated with the synthesis and secretion of fungal cellulases and xylanases, and a guide for P. oxalicum application in biotechnology.


INTRODUCTION
Soil microbiomes play an essential role in the cycles and balances of terrestrial ecosystems. A major population of organisms, filamentous fungi, break down organic matter to gain nutrients for survival, owing to their ability to synthesize and secrete plant-biomass degrading enzymes (PDEs; Gusakov, 2011). Natural enzyme production by fungi, however, is very low, restricting their industrial application. Although the isolation and characterization of PDEs has been performed, success has been limited, because of the diversity of PDEs, signal transduction specificity in response to carbon sources and the complex synergy of PDEs required to degrade plant biomass (Zhang and Zhang, 2013;Parameswaran et al., 2018;Schmoll, 2018;Antonieto et al., 2019;Srivastava et al., 2019). Breeding by genetic modification has been a potential strategy to improve fungal enzyme production. However, to best perform this strategy, it is necessary to understand the mechanisms responsible for PDE synthesis and secretion by fungi, specifically aspects associated with the expression of enzyme-encoding genes and their expressional regulation.
The effects of PDEs on plant-biomass degradation result from a co-expression network of genes mediated by carbon sources, including enzyme and other associated genes. Recently, the development of high-throughput sequencing technologies (Shendure et al., 2017) and their corresponding analysis tools, has resulted in decreasing sequencing costs ($0.012 per Mbp 1 ), thereby facilitating large-scale analyses of omics data that can uncover universal mechanisms at the molecular level. Weighted gene co-expression network analysis (WGCNA) is a systems biology method for describing correlation patterns among genes based on RNA-sequencing data and it is widely used for the construction of gene co-expression networks and identification of critical hub genes in several biological contexts, including yeasts, mammals and plants.
High expression of fungal PDE-encoding genes requires a corresponding inducer. At present, insoluble biomass [i.e., Avicel (AV), xylan and raw starch] and soluble sugars (i.e., D-xylose, cellobiose, lactose, and sophorose) are commonly used for inducing PDE production in filamentous fungi, including Trichoderma, Aspergillus, and Penicillium (Amore et al., 2013;Tani et al., 2014). Several studies postulated the inducer function of low molecular weight and soluble compounds from plant biomass. When cultivated on raw insoluble biomass, fungal cells initially produce a basal level of PDEs that hydrolyze biomass to generate a soluble inducer (Amore et al., 2013). Therefore, a potent inducer is essential for high PDE yields in filamentous fungi.
Transcription of fungal PDE-encoding genes is tightly controlled by specific transcription factors (TFs). Several TFs regulating the expression of PDE-encoding genes were identified, including the cellulase activator CLR-2/ClrB (Coradetti et al., 2012;Zhao et al., 2016), xylanase activator XlnR/Xyr1/XLR-1 (Amore et al., 2013;Li et al., 2015), amylase activator AmyR (Li et al., 2015;Benocci et al., 2017) and carbon catabolite repressor CreA/CRE1/CRE-1. These TFs can directly bind to the promoters of key PDE-encoding genes to perform their regulatory roles. Remarkably, XlnR negatively regulates cellulase production, while positively regulating xylanase production in P. oxalicum (Li et al., 2015). Both ClrB and XlnR dynamically control the transcriptional levels of major cellulase and/or xylanase genes under the induction of cellulose, but their relative importance depends on the species (Benocci et al., 2017). AmyR functions in the opposite way, in the expressional regulation between amylase and cellulase genes in P. oxalicum (Li et al., 2015). CreA/CRE1/CRE-1 negatively regulates the transcripts of almost all PDE-encoding genes and their regulatory genes (Li et al.,1 www.genome.gov/sequencingcostsdata 2015; Benocci et al., 2017). This suggests that the expression of different PDE-encoding genes is co-regulated, but shows diversity. In addition, several studies indicate that expression of PDE-encoding genes is also co-regulated with that of genes involved in asexual reproduction and secondary metabolites, such as BrlA and LaeA Li et al., 2016).
In this study, to find an efficient inducing carbon source for P. oxalicum, we investigated the effects of 23 carbon sources on cellulase and xylanase production in the P. oxalicum strain PoxKu70. Of these, two novel carbon sources [methylcellulose (MC) and 2-hydroxyethyl cellulose (HEC)] efficiently promoted enzyme production. Transcriptional profiling of the PoxKu70 was also performed, when cultured on AV, MC, HEC or WB. In addition, we employed WGCNA to analyze 24 time-course transcriptomes, with three biological replicates for each sample, from P. oxalicum under various culture conditions. The coexpression modules of cellulase and xylanase-encoding genes were identified and characterized in P. oxalicum. Moreover, the regulatory roles of modular TFs in the production of cellulase and xylanase were investigated.

Fungal Strains and Culture Conditions
The fungal strain used in this work is the highly efficient transformation strain PoxKu70 (#3.15650, China General Microbiological Culture Collection, CGMCC), derived from the P. oxalicum wild-type strain HP7-1 (CGMCC no. 10781) by deleting the gene PoxKu70 (Zhao et al., 2016). Fresh spores of PoxKu70 were statically incubated on potato-dextrose agar plates at 28 • C for 6 days, collected and resuspended in sterile water containing 0.1% Tween 80 to a final concentration of 1 × 10 8 per mL for further study.
For RNA-sequencing, the fresh spore suspension described above (500 µL) was inoculated into GLU medium (100 mL) (Yan et al., 2017) and placed in a shaker at 180 rpm for 24 h at 28 • C. The hyphae were filtered using a vacuum pump and washed once with sterile water (50 mL). Approximately 1 g of hyphae was transferred to basic medium (100 mL), separately supplemented with various carbon sources, including 1% (w/v) AV, WB, MC, HEV, and GLU, and incubated in a shaker at 180 rpm for 4-48 h at 28 • C. The medium with no carbon source (NC) and GLU were used as controls. Hyphae were harvested for total RNA extraction and the remaining medium (crude enzyme solution) was used for the measurements of enzymatic activities.

Total RNA Extraction
Total RNAs from P. oxalicum hyphae under different culture conditions were extracted using a TRIzol RNA Kit (Life Technologies, Carlsbad, CA, United States), as described previously (Yan et al., 2017).

RNA-Sequencing and Data Analysis
RNA-sequencing of P. oxalicum total RNAs was carried out on a BGISEQ-500 platform at BGI, Shenzhen, China. In brief, total RNAs were pretreated with DNase I to eliminate double-and single-stranded DNA contaminants, and the mRNA was purified after mRNA enrichment using oligo (dT), and attached to magnetic beads. Subsequently, the purified mRNA was fragmented and used as the template for the synthesis of first-strand cDNA with random hexamer-primed reverse transcription, followed by the synthesis of the second-strand cDNA. The generated cDNA was subjected to end repair and 5 -phosphorylation, as well as the addition of a single base A and ligation with a bubble adapter. The products were enriched by PCR amplification and purified. Double-strand DNA was denatured by heating and subsequently circularized to construct a library of single-strand circular DNAs for further sequencing.

Enzymatic Activity Measurement
Fungal cellulase and xylanase activities were measured as described previously (Yan et al., 2017).

Measurement of Intracellular Protein
The cultivated mycelium collected from 100 mL of culture was harvested, washed and ground into powder with liquid nitrogen. Protein extraction buffer [15 mL, containing NaCl 8.5, Na 2 HPO 4 2.2, NaH 2 PO 4 0.2, phenylmethylsulfonyl fluoride 0.87 and ethylenediaminetetraacetic acid 1.86 g/L, at pH 7.4] and 0.25 mm dia. glass balls (4 g) were added to the powder. The mixture was put on ice for 5 min, shaken for 1 min and put back on ice for 1 min, repeated in triplicate. The protein was 2 https://cran.r-project.org/web/packages/pheatmap/ separated via centrifugation at 13000 × g for 10 min. Protein concentration was determined using the Pierce TM detergentcompatible Bradford assay kit (Pierce Biotechnology, Rockford, IL, United States).

Measurement of Growth Curves of P. oxalicum
Hyphal weight from P. oxalicum was quantified as total intracellular protein content per 100 mL of culture. Fresh spores (10 8 ) of PoxKu70 were directly inoculated into medium (100 mL) containing MC, HEC, AV, WB, GLU, or NC, respectively. The inoculated media were cultivated at 28 • C for 24-72 h. The mixture of mycelia and/or insoluble carbon sources were collected from 100 mL of culture by vacuum filtration every 12 h, and then the entire mycelia collected were used for extraction of intracellular proteins as described above. The resulting intracellular protein content data were used for construction of growth curves for P. oxalicum.

WGCNA
The R package, WGCNA was used for WGCNA (Langfelder and Horvath, 2008). The average FPKM value of three biological replicates for each gene was calculated and used for further analyses. R package HCLUST 3 was used for detecting outliers of samples with the clustering method "average." The modules were identified using the dynamic tree cut method based on dissTOM hierarchical clustering and the other parameters, deepSplit = 2 and minModuleSize = 30. The modules having dissimilarity coefficients of 0.25 (cutHeight = 0.25) were merged to obtain final co-expression modules. The moduleTraitCor function was used to correlate (hemi-) cellulase activity traits with co-expression modules according to the tutorials of WGCNA (Langfelder and Horvath, 2016). Briefly, we used module Eigengenes from WGCNA package to process the expression profile in each module to generate the corresponding module eigengene, and subsequently related the module eigengenes to the activities of cellulases and xylanases of P. oxalicum strain PoxKu70 cultured on various carbon sources through calculating the correlation coefficients using moduleTraitCor. Significance analysis was carried out using corPvalureStudent. Heatmap showing the correlation coefficients and corresponding values was constructed using labeled heatmap. Values (pvalue ≤ 0.05 and correlation coefficients ≥ 0.4) indicate correlations between samples (Langfelder and Horvath, 2008). Relevance matrix calculations and significance level calculations of correlation coefficients between global genes were performed using the R-package Hmisc 4 .

Construction of Deletion Mutants
Deletion mutants were constructed in the P. oxalicum parental strain PoxKu70 according to the methods described by Yan et al. (2017). The primers used for mutant construction are listed in Supplementary Table S1.

Availability of Data and Materials
All transcriptomic data are available in the GEO database with accession number GSE133258 on NCBI. The DNA sequences for POX01118, POX01474, and POX01678 have been deposited in the GenBank database (accession numbers MN529555-MN529557), respectively.

RESULTS AND DISCUSSION
Enzyme Activities From P. oxalicum

Induced by Various Carbon Sources
To screen for strong inducers of cellulase activity in the PoxKu70, 23 carbon sources were chosen (Supplementary Table S2). NC was used as a control. Activities of four cellulases [filter paper cellulase (FPase), carboxymethylcellulase (CMCase), p-nitrophenyl-β-cellobiosidase (pNPCase) and p-nitrophenyl-βglucopyranosidase (pNPGase)], plus xylanase were measured after 2-4 days of culture after transfer from glucose. The activities of cellulases and xylanases showed diversity under the induction of different carbon sources. For example, FPase activities of PoxKu70 with more than 0.01 U/mL were detected when cultured on AV, CMC, HEC, MC, WB, α-cellulose and Sigmacell cellulose (SC50; Sigma-Aldrich, Darmstadt, Germany) for 2 days after a transfer from GLU, or on additional lactose and cellobiose (Sigma-Aldrich) for 4 days. Of these, FPase activity on MC was the highest, at 0.16 U/mL at 2 days and 0.36 U/mL at 4 days, followed by WB (0.08 U/mL at 2 days) and α-cellulose (0.20 U/mL at 4 days), respectively. CMCase activities of the PoxKu70 ranged from 0.02 to 4.31 U/mL when cultivated on lactose, CMC, cellobiose, HEC, MC, α-cellulose, AV, WB, D/Larabinose, or Sigmacell cellulose, but they could not be detected on other carbon sources. PoxKu70 cultivated on MC possessed the highest activity of CMCase (i.e., 1.08 U/mL at 2 days and 4.31 U/mL at 4 days), followed by AV (2.91 U/mL at 4 days) and WB (1.75 U/mL at 4 days), respectively. Both pNPCase and pNPGase activities were low, ranging from 0.01 to 0.89 U/mL on lactose, CMC, cellobiose, HEC, MC, xylose, α-cellulose, AV, WB and Sigmacell cellulose. WB was the most efficient carbon source to induce pNPCase and pNPGase activities of P. oxalicum, among the tested carbon sources. Xylanase activity in PoxKu70 was 0.02-16.07 U/mL on 12 carbon sources, including lactose, CMC, cellobiose, HEC, MC, xylose, α-cellulose, AV, WB, D/Larabinose and Sigmacell cellulose, except for D-arabinose at 2 days. Xylanase activity on α-cellulose was higher than that of other carbon sources. PoxKu70 hardly produced any cellulases or xylanases when cultured on NC. Overall, the inducibility of FPase and CMCase production by MC was strongest at 2-4 days. WB was the most efficient carbon source to induce pNPCase and pNPGase production and α-cellulose for xylanase production, among the tested carbon sources (p < 0.05; Figure 1).
Methylcellulose and HEC are novel carbon sources for induction of cellulase and xylanase production from P. oxalicum. Compared with the already known WB and AV, induction by MC was stronger, specifically for FPase and CMCase, while slightly lower for HEC. MC and HEC are derived from cellulose through methylation and hydroxyethylation, respectively, and are widely used in tissue engineering (Dutta et al., 2019). MC mainly induced the secretion of EGs that randomly cleave internal β-1,4-glycosidic bonds of cellulose chains to generate numbers of chain ends for CBHs, thereby resulting in increased FPase activity. By contrast, HEC stimulated all types of cellulases, including CBHs, EGs and BGLs, and increased FPase activity, dependent on a strong synergistic interaction. Remarkably, strong inducers for Tricherdoma and Aspergillus, such as lactose (Amore et al., 2013), showed low induction of cellulase and xylanase production in P. oxalicum, probably through repressing the expression of cellulase and xylanase-encoding genes, but this needs to be confirmed. These data also imply that different signal transduction pathways and regulation mechanisms for cellulase and xylanase-encoding genes exist in P. oxalicum.

Effects on P. oxalicum Growth by the Five Chosen Carbon Sources Above
Five carbon sources including two novel inducers MC and HEC, commonly used inducers AV, WB and repressor GLU, were selected for further study. The growth curves of PoxKu70 were determined after direct inoculation into medium containing AV, WB, MC, HEC, and GLU, respectively, with an NC control. Here the intracellular protein contents measured from the entire mycelia collected from the inoculated media were used for construction of the growth curves. Growth of PoxKu70 cultured on GLU, WB, and AV, was faster than that on HEC, MC, and NC, except on AV before 24 h (Figure 2). Biomass accumulation of PoxKu70 on AV at 24 h was the same as that on HEC, NC, and MC. During the early period of culture (before 36 h), PoxKu70 grew similarly on WB and on GLU since WB contains some free GLU, but growth on WB was slower after 36 h, since free GLU would have been consumed and GLU then needed to be released via WB hydrolysis by cellulases or amylases. On AV, PoxKu70 grew slower than on WB before 36 h, faster after 64 h, but slower than on GLU during the whole culture time. The growth of PoxKu70 was similar, when cultured on HEC, MC or NC from 24 to 36 h, while decreased on MC and increased on HEC compared with on NC after 36 h (Figure 2).
WB is rich in polysaccharides (i.e., starch and cellooligosaccharides), simple sugars such as GLU and crude proteins (Sun et al., 2008;Sardari et al., 2019). When cultured on WB, P. oxalicum initially grew on the GLU in WB, as well as it did on GLU medium. Once the GLU was exhausted, the P. oxalicum cells would have been forced to produce cellulases and/or amylases to degrade cello-oligosaccharides and/or starch into GLU, resulting in a growth delay (Sun et al., 2008). Avicel is a crystalline cellulose that cannot be directly utilized by fungal cells (Yang et al., 2011), but can be digested into GLU by cellulases from P. oxalicum. Therefore, mycelium accumulation gradually increased during the culture time, dependent on GLU release. During the early culture period, the growth of P. oxalicum on HEC was similar to that on NC because the inoculated hyphae had only stored energy available, but sped up after nutrients were released from HEC degradation. Surprisingly, MC inhibited fungal growth, which FIGURE 1 | Activities of crude enzymes produced by P. oxalicum strain PoxKu70 cultured on various carbon sources. Crude enzymes were collected from the PoxKu70 cultivated for 2-4 days after a transfer from GLU. * , p ≤ 0.05; * * , p ≤ 0.01 indicate the significance of differences between the results on various carbon sources and WB. CMC, carboxymethylcellulose; HEC, 2-hydroxyethyl cellulose; MC, methyl cellulose; WB, wheat bran; AV, Avicel; D-arab, D-arabinose; L-arab, L-arabinose; SC50, Sigmacell cellulose; NC, no carbon source; GLU, glucose; FPase, filter paper cellulase; CMCase, carboxymethylcellulase; pNPCase, p-nitrophenyl-β-cellobiosidase; pNPGase, p-nitrophenyl-β-glucopyranosidase.
Frontiers in Microbiology | www.frontiersin.org PoxKu70 were inoculated into medium containing each carbon source and cultivated for 24 to 72 h. Total intracellular proteins from all the P. oxalicum hyphae in the culture were extracted and measured. The lowercase letters indicate significant difference at p < 0.05 probability level, assessed by Student's t-test. HEC, 2-hydroxyethyl cellulose; MC, methyl cellulose; WB, wheat bran; AV, Avicel; NC, no carbon source; GLU, glucose.
might be attributable to toxicity of the released methyl-glucose or oligosaccharides from MC degradation.

Carbon Source-Specific Transcriptional Patterns of P. oxalicum
To further investigate carbon source-specific transcriptional patterns of P. oxalicum growth on the five carbon sources WB, HEC, AV, MC, and GLU, time-course transcriptomes from PoxKu70, cultured at four time points (4, 12, 24, and 48 h), after a transfer from GLU, were RNA-sequenced, with NC as control. Each sample had three biological replications. A total of 72 transcriptomes was sequenced. A total of 1,548,422,794 clean reads was generated, with an average length of 50 bp per read (Supplementary Table S3), and were mapped to the genome of P. oxalicum wild-type strain HP7-1 (Zhao et al., 2016), resulting in an average 98.2% mapping rate. Correlation coefficients were high (R > 0.82) for the three biological replicates produced under each set of culture conditions ( Supplementary  Figures S1-S4). Furthermore, a principal components analysis revealed no outliers among the three biological replicates for each sample (Supplementary Figure S5), indicating that all the transcriptomes were reliable for further analyses.
We comparatively analyzed 72 time-course transcriptomes described above using DESeq software (Love et al., 2014). Using p-value ≤ 0.05 and | log2 fold change| ≥ 1 as thresholds, differentially expressed genes (DEGs) were screened in transcriptomes of P. oxalicum on five carbon sources (AV, WB, HEC, MC, and GLU) compared with on NC as the control, respectively. Each comparison ( PoxKu70_AV vs.
Carbon sources are able to induce the genes encoding proteins involved in carbon metabolism and the tricarboxylic acid (TCA) cycle (He et al., 2013). Here, comparative transcriptomic analyses indicated that 60 DEGs were annotated to be involved in glycolysis/gluconeogenesis and the TCA cycle in P. oxalicum. These DEGs appeared at least once, in all comparisons between the transcriptomes expressed during growth on the five carbon sources and NC. These proteins encoded by DEGs included the major rate-limiting enzymes, hexokinase POX01079, POX04253, POX04475, and POX05787; phosphofructokinase POX08888; pyruvate kinase POX02744 in glycolysis, and citrate synthases POX01075 and POX04356; isocitrate dehydrogenases POX06270 and POX09065; oxoglutarate dehydrogenases POX00538 and POX09663 in the TCA cycle (Supplementary Figure S10). Of the 60 DEGs, 7-20, 10-34, 11-15, 1-7, and 12-38 DEGs were susceptible to induction by growth on AV, WB, MC, HEC, and GLU for 4-48 h, to various degrees, respectively, compared with NC ( Figure 5). These results suggested that the expression of glycolysis/gluconeogenesis-and TCA cycle-related genes are most strongly induced by GLU and WB, followed by MC and AV. Most of the DEGs involved in glycolysis/gluconeogenesis and the TCA cycle were up-regulated in response to AV, MC, WB, and HEC, especially 20 of 34 DEGs with WB at 4 h compared with NC, but generally, much less up-regulated than in response to GLU (Figure 5). Most of 60 DEGs showed active alteration during the early period (up to 12 h) of growth on WB and GLU. Surprisingly, the transcriptional levels of those DEGs encoding rate-limiting enzymes showed variance under different conditions, for example, genes POX01079 and POX05787 were down-regulated on AV, WB, MC, and HEC compared with FIGURE 3 | Venn diagram indicating numbers of unique and shared DEGs in transcriptomes from P. oxalicum strain PoxKu70 in the presence of different carbon sources (AV, MC, HEC, WB, and GLU) compared with NC. DEGs were selected using the thresholds of p-value ≤ 0.05 and | log2 fold change| ≥ 1. DEGs, differentially expressed genes; AV, Avicel; MC, methyl cellulose; HEC, 2-hydroxyethyl cellulose; GLU, glucose; NC, no carbon source.
NC, but up-regulated on GLU at 4 h. At 12 h, however, their transcriptional levels showed no alteration (Figure 5).
The detected DEGs included other genes possibly involved in synthesis of cellulase and xylanase, such as TFs, transporters and G protein-coupled receptor (GPCR) genes. Comparative transcriptomic analyses showed that 283 DEGs encoded putative TFs. Approximately 70% of these detected TFs belong to the zinc finger family, followed by winged helix repressors (8.4%). Their transcriptional levels varied in response to the different carbon sources and the duration of induction, for instance, 44 and 165 DEGs on MC and WB at 4 h compared with that on NC, respectively, and 34 and 51 at 12 h. Screening for the top ten upregulated/downregulated TF genes for each comparison detected a total of 163, including several known regulatory genes, i.e., POX01167/PoxCxrA, POX04420/PoxCxrB, POX01960/PoxClrB, POX02484, POX03890/PoxAmyR, POX08292/PoxMBF1, and POX06534/BrlA Li et al., 2015;Zhao et al., 2016Zhao et al., , 2019Yan et al., 2017; Figure 6).
A total of 271 DEGs encoding proteins were annotated as ABC transporter-like (IPR003439), sugar/inositol transporter (IPR003663), or major facilitator superfamily (IPR011701). Numbers of transporter DEGs in response to AV and MC increased, after longer induction duration compared with NC, but numbers decreased in response to GLU and WB, and hardly changed in response to HEC. Homologous analysis with known GPCR in filamentous fungi (Affeldt et al., 2014;Martín et al., 2019) detected 18 DEGs encoding GPCRs. Their expression was active under induction by WB and GLU during the early period (Figure 6). The comparative transcriptomic analysis detected many more DEGs, but most of them were not related to cellulase and xylanase production in filamentous fungi; this needs to be studied further.

Identification of Co-expression Modules of Cellulase-and Xylanase-Encoding Genes in P. oxalicum
We hypothesized that the genes mediated by carbon sources would form tight modules, also called associations, with each other, in gene co-expression networks. The 72 transcriptomes obtained were subjected to WGCNA (Langfelder and Horvath, 2008). Sample outlier detection was evaluated first. The average values of FPKMs generated by three biological repeats for each sample were used as input data for a WGCNA. Sample outlier detection revealed that MC48, transcriptome from P. oxalicum cultured on MC for 48 h, had the highest outlier potential. The outlier potential of MC48 could be attributed to the efficient induction of gene expression by MC at 48 h, because the enzyme extracts of P. oxalicum showed the highest cellulase and xylanase activities, except for pNPGase (Supplementary Figure S11). Therefore, MC48 needed to be retained for further WGCNA.
Subsequently, the pickSoftThreshold function was used to calculate the soft thresholds. When soft-power was set at 24, the Scale Free Topology Model Fit R 2 was ≥ 0.9 (Figures 7A,B). To confirm that the topology of the gene co-expression network was scale-free, soft connectivity was determined. The correlation coefficient between the log (k) of the node number with connectivity K and the log [p (k)] of the probability of node occurrence was ≥ 0.9 (Figure 7C), suggesting that the topology of the network trended scale-free when the soft-power was 24.
A total of 45 modules were identified using the Dynamic Tree Cut method Horvath, 2007, 2008; Figure 7D). Further merging of modules generated 17 modules with set module dissimilarity coefficients of ≤ 0.25 ( Figure 7D). Members of these 17 modules showed variances, ranging from 36 to 2,623 (Figure 8). Of these, MEgrey contained genes that were not co-expressed. The genome of P. oxalicum contained 35 genes annotated as cellulases and XYNs (Zhao et al., 2016) that were further analyzed. They were divided into seven modules, MEivory, MEmidnightblue, MEorange, MEroyalblue, MEturquoise, MEgrey and MElightcyan1. WGCNA analysis of T. reesei grown on sugarcane bagasse identified 28 highly connected gene modules and one of these modules contained the most representative core of cellulolytic enzymes, with their regulators and transporters (Borin et al., 2018).
Among the nine putative TF genes, four genes, POX01960/ClrB, POX02071/ClrB-2, POX03837/ClrA and POX05324/XlnR, encoded Zn2Cys6 proteins that are involved in the regulation of cellulase and/or xylanase production in filamentous fungi (Coradetti et al., 2012;Li et al., 2015;Zhao et al., 2016). ClrB up-regulates cellulase and xylanase production, except for pNPGase (Zhao et al., 2016), while a deletion of ClrB-2 partially reduces the cellulase activity of P. oxalicum (Li et al., 2015). XlnR down-regulates cellulase production while up-regulating xylanase production (Li et al., 2015). Both ClrB and XlnR dynamically control the transcriptional levels of major cellulase and/or xylanase genes under the induction of cellulose, but their precise roles differ depending on the species concerned , the soft-powder was selected as 24, resulting in a scale-free gene co-expression network topology. (D) Identification of functional modules illustrated with various colors. Cluster dendrograms were generated using hierarchical clustering. Original modules were identified using the Dynamic Tree Cutting method, and merged modules were generated according to the correlations of modules, with module dissimilarity coefficients of ≤ 0.25. (Benocci et al., 2017). ClrA/Clr-1 is not only required for the expression of major cellulase and xylanase-encoding genes, but is also necessary for the expression of ClrB/Clr-2 (Coradetti et al., 2012). However, whether the five remaining TF genes POX01118, POX01474, POX01678, POX04007 and POX07871 participate in the regulation of cellulase and xylanase production of P. oxalicum, needs to be further studied.
Among the 15 transmembrane proteins, POX06051, encoding a key cellodextrin transporter POX06051/CdtC, is required for cellobiose consumption and the expression of major cellulase genes in P. oxalicum, in combination with another transporter, POX05915/CdtD (Li et al., 2013), suggesting the crucial role of sugar transporters in the induction of cellulase gene expression by insoluble cellulose. In addition, the inductive roles of sugar transporters exhibit redundancy (Li et al., 2013). Here, seven of the remained ten genes encoding transmembrane proteins were annotated as sugar/inositol transporters (IPR003663) and/or members of the major facilitator superfamily (IPR011701), without well-characterized biological functions. POX00907, POX01263, and POX07722 were predicted to encode a cytochrome P450, an integral membrane protein and a bicarbonate transporter, respectively. Two candidate classic G protein-coupled receptors POX01263 and POX03930 containing seven transmembrane domains (TMDs) had a specific feature of classic GPCRs (Martín et al., 2019). However, POX03930 was annotated as a RAT1-like protein that is involved in fungal detoxification (Soustre et al., 1996).
The secondary module MEroyalblue contained 157 members (Supplementary Table S5), including seven cellulase and XYN genes (a bgl POX08882, an eg POX07535/Cel12A and five xyns POX04274, POX05916, POX06601, POX06783/Xyn11A and POX08990) and 22 other CWDE genes. Notably, most of the xylanolytic genes in the HP7-1 genome were included in MEroyalblue, such as five xyns described above and eight β-xylosidase genes (POX00007, POX01914, POX01921, POX05540, POX06599, POX06600, POX07441, and POX07891). Only two TF genes POX03888/PrtT and POX06865 were included. POX03888/PrtT, a transcription activator of extracellular proteases in filamentous fungi, could directly regulate the expression of many peptidase genes and several amylase genes including α-amylase, glucoamylase and α-glucosidase genes. POX03888/PrtT also indirectly regulated the expression of cellulase genes, which may be related to nutrient limitations (Chen et al., 2014). However, to our knowledge, there is no report about its regulation of xylanolytic genes in P. oxalicum. MEroyalblue also contained 19 transmembrane proteins with 2-14 TMDs, most of which were annotated as sugar/inositol transporter (IPR003663) and/or major facilitator superfamily (IPR011701), including cellodextrin transporter POX05915/CdtD. The remaining 92 of 157 members in module MEroyalblue included 38 genes involved in metabolism (i.e., amino acid metabolism, carbohydrate metabolism, energy metabolism, lipid metabolism and enzyme families) and 54 unknown genes through KEGG annotation. Notably, four genes POX02832, POX02775, POX05829, and POX07164, encoded putative aldose 1-epimerase GalM, alcohol dehydrogenase, alcohol dehydrogenase and triosephosphate isomerase TpiA, respectively, which participate glycolysis/gluconeogenesis. However, these annotated genes were not previously known to be associated with cellulase and/or xylanase production in P. oxalicum.
The remaining five modules MEmidnightblue, MEorange, MEturquoise, MEgrey, and MElightcyan1 contained a few unimportant cellulase and XYN genes, and their members accounted for approximately 80% of all genes annotated in the genome of HP7-1 (Zhao et al., 2016), which should not therefore, be co-expressed with the major cellulase and XYN genes.
In order to understand the significance of the identified modules toward cellulase and xylanase production of P. oxalicum, their correlation coefficients were calculated. Using correlation coefficient > 0.4 and p-value < 0.05 as thresholds, the relevant modules to cellulase and xylanase production of P. oxalicum were screened and characterized. The correlation strengthens with the increase in the correlation coefficient. Module MEivory was most relevant to filterpaper cellulase and CMCase, with correlation coefficients of 0.92 (p-value 3e−10) and 0.76 (p-value 2e−5), respectively, followed by MEdarkgrey (0.43, p-value 0.04, and 0.44, p-value 0.03, respectively) and MElightsteelblue1 (0.47, p-value 0.02, and 0.51, p-value 0.01, respectively). MElightcyan1 showed low negative correlations (−0.43, p-value 0.04 and −0.47, p-value 0.02) with FPase and CMCase production, respectively, whereas no module was found to be correlated with pNPCase or pNPGase production. Only MEivory correlated with xylanase production of P. oxalicum (0.72, p-value 3e−10; Figure 9A). These results were confirmed through gene and module significance analyses (Langfelder and Horvath, 2008; Figures 9B,C).
In addition, a relationship assay between module eigengenes indicated that MEivory had relatively tight correlations with MElightsteelblue1, MElightcyan1, and MEdarkgrey, compared with other modules, having correlation coefficients of 0.51, 0.47, and 0.41 (Figure 10)  (Ruger-Herreros and Corrochano, 2020). For example, GLU starvation induces hydrolase production and maturation, while GLU richness enhances mycelial biomass and conidiation. The secretion of extracellular hydrolases supports conidia formation by releasing nutrients (Emri et al., 2018) such as GLU, resulting in the co-regulation of fungal conidiation and hydrolase production. Two key members of the central regulatory pathway, BrlA and AbaA, tightly control fungal sporulation, acting in concert with fluffy low BrlAs (FLBs), such as FlbB and FlbD (Park and Yu, 2012). BrlA down-regulated cellulase activity of P. oxalicum , while AbaA acted positively (data not shown). Proteins AygA, AbrA, Alb1/wA, ArpA, and ArpB are involved in pigment biosynthesis (Mayorga and Timberlake, 1990;Park and Yu, 2012), and NimX, acyclin-dependent kinase, controls fungal cell division (McGuire et al., 2000). TmpA is a putative membrane flavoprotein that regulates asexual development in Aspergillus nidulans (Soid-Raggi et al., 2006). Whether these proteins are involved in cellulase and xylanase production of P. oxalicum needs to be studied further.
In MElightcyan1, the gene POX07254 encodes an essential TF, CreA, that mediates fungal carbon catabolite repression in the presence of the favored carbon source, GLU. CreA inhibits the expression of almost all the cellulase and xylanase-encoding genes and their regulatory genes, including ClrB and XlnR, induced by insoluble cellulose (Antonieto et al., 2014;Li et al., 2015;Huberman et al., 2016). In contrast, the other genes in this module have neither been found to be associated with the regulation of cellulase and xylanase-encoding genes in fungi, nor as genes in module MElightsteelblue1.

Identification of Highly Connected (Hub) Genes in MEivory
The co-expression network of genes in MEivory was constructed using Cytoscape 3.6.1 software (Shannon et al., 2003), and the hub genes were investigated in the constructed network. The 20 largest hub genes in MEivory were POX03368, POX06079, POX06051/CdtC, POX03367, POX03641, POX02071, POX01166/Cel5B, POX05571/Cel7B, POX05570/Cel45A, POX01896/Cel5C, POX01646, POX02739, POX01961, POX01960/ClrB, POX01833, POX03062, POX05587/Cel7A-2, POX08897/AA9, POX07524, and POX07415, having 34 to 48 coupled genes (Figure 11). Most of them encoded proteins that were CAZymes, including the essential cellulases (CBH POX05587/Cel7A-2; EGs POX01166/Cel5B, POX05571/Cel7B, POX05570/Cel45A and POX01896/Cel5C; BGLs POX01646, POX03062, POX06079, and POX03641; Auxiliary activity POX08897/AA9), and their regulators such as POX01960/ClrB and POX02071/ClrB-2. It is not surprising that these cellulase genes and their regulatory genes in filamentous fungi including P. oxalicum, T. reesei, and Neurospora crassa, are clustered together because they are co-induced in response to most plant biomass substrates such as WB, sugarcane bagasse and Avicel (Craig et al., 2015;Yan et al., 2017;Borin et al., 2018). Under fungal growth conditions, in the presence of cellulose, cellulose is first digested into cellodextrins by constitutive levels of extracellular CBHs and EGs. The generated cellodextrins are transported into the fungal cells by the cellodextrin transporter CdtC, then either degraded into GLU by intracellular β-glucosidases, such as POX06079, or act as an inducer of cellulase and xylanase gene expression (Li et al., 2015). The proteins POX01646, POX01833, POX01961, POX02739 and POX07524, were annotated as β-xylosidase, α-mannosyltransferase, β-mannosidase, carbohydrate acetylesterase and α-L-rhamnosidase, respectively, and may assist the major cellulases with degradation of plant cell walls. POX03368 shared a 33.72% identity with MreA from Aspergillus oryzae (accession number BAB13480.1) that encodes an isoamyl alcohol oxidase. Isoamyl alcohol oxidase catalyzes the formation of isovaleraldehyde as the main component of mureka, producing the off-flavor (Yamashita et al., 2000). POX03367 was annotated as an ATP-dependent Clp protease, ATP-binding subunit ClpX. ClpX is a protein chaperone of ClpP proteases, which are responsible for protein homeostasis through the degradation of damaged or unfolded proteins, as well as for the conditional degradation of functional proteins in response to environmental signals (Stahlhut et al., 2017). The function of ClpXP has been studied broadly in prokaryotes but rarely reported in fungi. Whether POX03367 and POX03368 are involved in the regulation of cellulase and xylanase-encoding gene FIGURE 11 | Co-expression network of the highly connected genes in module MEivory. The co-expression network of genes in MEivory was constructed using Cytoscape 3.6.1 software. This network only shows the genes having correlation coefficients greater than 0.4. Cycle size indicates the number of connected genes. Green, CWDEs; Blue, the remaining CAZymes; Light-green, TFs; Dark-orange, proteins involved in protein homeostasis; Orange, other annotated proteins; Pink, putative transporters; Light-blue, hypothetical proteins; Purple, transferases. expression levels in fungi remains unknown, and needs to be further studied.

Kinetic Analysis of the Expression Levels of MEivory Members
The overall expression levels of these 120 genes were induced by complex carbon sources (i.e., AV, MC, WB, and HEC), but the levels varied depending on the carbon source. Conversely, their transcriptional levels in both the presence of GLU and NC were low or inhibited ( Figure 12A). Comparative analyses revealed that the expression levels of MEivory members increased between 4 and 24 h of induction by AV and then decreased, while remaining increased over time by MC. In contrast, their expression only increased from 4 to 12 h and then stabilized in the presence of HEC or WB. The expressional levels of MEivory in the presence of MC and AV were higher than those on HEC and WB, especially at later stages of culture (Figure 13). The overall expression changes over time should be attributed to the distinct structure and chemical composition of different biomass substrates (Wang et al., 2015). AV has a high crystallinity, in which large macrofibrils are aggregated together, resulting in a tightly packed structure, which limits the penetration of cellulolytic enzymes. Therefore, degradation of AV is based on a layer-by-layer model, with synergistic interactions between CBHs, EGs and BGLs (Yang et al., 2011;Bornscheuer et al., 2014). During the early stage of AV induction, P. oxalicum stimulates the expression of genes in MEivory to degrade AV and promote GLU uptake, which increases over time. However, during the later stages, the secreted cellulolytic enzymes were enough for the degradation of the remaining substrate, or fungal cells might  have required little GLU as a result of maturation or death, which contributed to the reduction in gene expression. MC and HEC are derived from cellulose, through random methylation and hydroxyethylation, respectively. Theoretically, their digestion also requires synergistic action by several cellulolytic enzymes. However, the methylated and hydroxyethylated groups should inhibit the binding of enzymes to these substrates, which might lead to a different enzyme composition. In addition, methylated and hydroxyethylated cello-or mono-oligosaccharides might enter into fungal cells to differentially induce gene expression, but this needs confirmation. In contrast, WB is composed predominantly of starch, non-starch polysaccharides such as cello-oligosaccharides and GLU, and crude proteins (Sardari et al., 2019). Cello-oligosaccharides are required for cellulase and xylanase production in P. oxalicum (Sun et al., 2008). However, how WB induces the expression of genes encoding cellulolytic enzymes has been unclear up to now. GLU is known to repress the secretion of cellulolytic enzymes in filamentous fungi through promoting carbon catabolite repression (CCR). In the presence of GLU, CCR is regulated mainly by transcription factor CreA, that represses the expression of genes encoding plant cell wall-degrading enzymes and their regulatory genes. The regulation of CreA partially depends on its intracellular localization and phosphorylation (Huberman et al., 2016). Under growth induction conditions, CreA forms a CreA-SsnF-RcoA repressor complex that may promote CreA degradation, connected to an SCF ubiquitin ligase complex, via GskA protein kinase, whereas, when supplied with GLU, CreA disconnects from the ubiquitin ligase complex and is then transported into the nucleus of A. nidulans (de Assis et al., 2018).
The transcripts of two key TF genes, POX01960/ClrB and POX05324/XlnR, directly regulate the expression of major cellulase and xylanase-encoding genes in P. oxalicum. Transcription of POX01960/ClrB was higher than that of POX05324/XlnR during the whole culture time on each carbon source (Figure 12B), suggesting that POX01960/ClrB plays a much greater role in the regulation of cellulase and xylanaseencoding gene expression in P. oxalicum than POX05324/XlnR. This finding is in contrast to Xyr1 and XlnR in T. reesei (Mach-Aigner et al., 2008;Borin et al., 2018), which was confirmed through the previous genetic analysis (Li et al., 2015;Zhao et al., 2016).
In addition, all 12 DEGs increased their transcriptional levels to various degrees under the induction of AV, MC, HEC and WB, with log2 fold changes from 1.1 to 10.3, except POX01833 on WB at 4 h. Conversely, most of them showed a significant reduction on GLU (Figure 12A).
In addition to these 12 DEGs that were induced commonly by biomass, the expression levels of several specific DEGs showed significant alterations that were dependent on the duration of induction and/or different substrates. For example, the expression levels of two key cbh genes, POX05587/Cel7A-2 and POX04786/Cel6A, were activated at 4 h on AV, HEC and MC, but not on WB, compared with NC. The regulatory gene POX01960/ClrB was induced after 24 h on AV and MC, as well as at 4 and 48 h on HEC, but was inhibited at 4 h on WB. The transcription levels of POX05324/XlnR on all five carbon sources showed no significant differences from NC, except at 24 h on MC (Figures 12, 14).
Consequently, the expression levels of major cellulase and xylanase genes in P. oxalicum were simultaneously induced, but at different expression levels, thereby resulting in different enzyme production levels. The expression changes may be attributed to different sites of action of the different cellulolytic enzymes and the distinct chemical composition of different biomass substrates (Wang et al., 2015). EGs randomly attack internal β-1,4-glycosidic bonds in cellulose chains to generate numbers of cello-oligosaccharides with exposed ends. CBHs release cellobiose from both ends of cello-oligosaccharides. Finally, BGLs break cellobiose and cellooligosaccharides into GLU (Wang et al., 2013). Atomic force microcopy indicated that the rough surface of the original AV became smoother after enzymatic hydrolysis, which results from the expression of the essential cbhs and egs during the early induction period (Yang et al., 2011;Bornscheuer et al., 2014); the AV derivatives HEC and MC behaved similarly. In contrast, WB is composed predominantly of starch, nonstarch polysaccharides and crude proteins. The non-starch polysaccharides included cello-oligosaccharides, arabinoxylans and β-(1,3) (1,4)-glucan (Sun et al., 2008). Therefore, when cultivated on WB, P. oxalicum could use GLU from the degradation of other polysaccharides, such as starch, during the early induction period.

Identification of the Uncharacterized TFs for Regulating Fungal Cellulase and Xylanase Activities
Weighted gene co-expression network analysis revealed that the expression of nine TF genes was co-expressed with major cellulase and xylanase-encoding genes in P. oxalicum, of which four were known. We further constructed the deletion mutants of three genes, POX01118, POX01474, POX01678 (Supplementary Figure S13) and investigated their regulatory roles in cellulase and xylanase production. Unfortunately, deletion of POX04007 and POX07871 that encoded C2H2and Zn2Cys6-zinc finger proteins failed. Enzymatic activity assays indicated that mutants POX01118 and POX01474 lost approximately 30% of their cellulase and xylanase activities under some induction conditions, compared with PoxKu70. For example, POX01118 had decreased FPase activity by 31.6-36.8% after 4 day of Avicel induction, while POX01474 had decreased xylanase activity by 22.0-29.8% at 4 day. Conversely, POX01678 had increased cellulase and xylanase activities (by approximately 30%) corresponding to FPase, CMCase, pNPCase and pNPGase, as well as xylanase activity (p < 0.05; Figure 15).

CONCLUSION
In this study, we identified a novel inducer (MC) for fungal cellulase and xylanase production and further comparatively analyzed the transcriptional patterns of P. oxalicum in response to MC, HEC, WB, AV, GLU, and NC. Moreover, through WGCNA analysis, the co-expression module MEivory that contains the cellulase and xylanase-encoding genes of P. oxalicum was identified. The transcriptional levels of MEivory members were dependent on the inducing carbon source and the induction time. Finally, three novel TFs POX01118, POX01474, and POX01678 were found to regulate the cellulase and xylanase production of P. oxalicum. These findings will facilitate understanding of the mechanisms required for the synthesis and secretion of fungal cellulases and xylanases and will provide a guide for P. oxalicum application in biotechnology.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the GEO database with accession number GSE133258 on NCBI. The DNA sequences for POX01118, POX01474 and POX01678 have been deposited in the GenBank database (accession numbers MN529555-MN529557).

AUTHOR CONTRIBUTIONS
J-XF conceived, supervised this study, and revised the manuscript. SZ codesigned and co-supervised this study, and revised manuscript. C-XL carried out bioinformatic analyses, mutant construction and enzymatic activity assay, and manuscript preparation. X-ML was involved in preparation of experimental materials. All authors read and approved the final version of the manuscript.