Construction of Light-Responsive Gene Regulatory Network for Growth, Development and Secondary Metabolite Production in Cordyceps militaris

Simple Summary Cordyceps militaris is an edible fungus that has been long used in traditional medicine. A large body of research has provided evidence of the medicinal properties of C. militaris extract and its demand has been increasing over the years. This study aims to construct and understand the light-responsive gene regulatory network of the fungi by combining the transcription factor (TF)-target gene interactions with the transcriptomic analysis of C. militaris under a light-programming condition. The study identified several key TFs and their gene targets that regulate growth, development and secondary metabolite production in the fungi under specific light conditions. Abstract Cordyceps militaris is an edible fungus that produces many beneficial compounds, including cordycepin and carotenoid. In many fungi, growth, development and secondary metabolite production are controlled by crosstalk between light-signaling pathways and other regulatory cascades. However, little is known about the gene regulation upon light exposure in C. militaris. This study aims to construct a gene regulatory network (GRN) that responds to light in C. militaris. First, a genome-scale GRN was built based on transcription factor (TF)-target gene interactions predicted from the Regulatory Sequence Analysis Tools (RSAT). Then, a light-responsive GRN was extracted by integrating the transcriptomic data onto the genome-scale GRN. The light-responsive network contains 2689 genes and 6837 interactions. From the network, five TFs, Snf21 (CCM_04586), an AT-hook DNA-binding motif TF (CCM_08536), a homeobox TF (CCM_07504), a forkhead box protein L2 (CCM_02646) and a heat shock factor Hsf1 (CCM_05142), were identified as key regulators that co-regulate a large group of growth and developmental genes. The identified regulatory network and expression profiles from our analysis suggested how light may induce the growth and development of C. militaris into a sexual cycle. The light-mediated regulation also couples fungal development with cordycepin and carotenoid production. This study leads to an enhanced understanding of the light-responsive regulation of growth, development and secondary metabolite production in the fungi.

called protein-binding microarray (PBM), to identify TF binding sequences. Weirauch et al. [29] analyzed data from the PBM assays of over 1000 TFs to determine DNA binding specificity represented by position frequency matrices (PFMs). The PFMs were then used to infer PFMs of other TFs from diverse eukaryotes, including C. militaris, based on the similarity of protein sequences [29]. To construct a gene regulatory network (GRN), the Regulatory Sequence Analysis Tools (RSAT) allows users to scan the upstream sequences of genes with the PFMs using Markov chain-based background models to identify TF binding domains [30]. The predicted TF-target binding domain interactions from RSAT are then used to infer the GRN.
This study aims to construct the GRN in C. militaris that responds to light. First, the collected PFMs of TFs of C. militaris were used to predict all possible TF binding sites in the promoter region of all genes in the C. militaris genome. Next, the interaction profiles between TFs and target genes were used to create a genome-scale GRN, which can be visualized as a network. Then, differentially expressed genes (DEGs) under light and dark conditions from transcriptome analysis were integrated into the network to extract a lightresponsive sub-network. The light-responsive GRN revealed the transcriptional regulation of many developmental and carotenoid and cordycepin biosynthesis genes by five key TFs, enabling more understanding of the regulation of light-controlled development and secondary metabolite production in C. militaris.

Materials and Methods
The overall workflow of the present study is shown in Figure 1 and is detailed as follows.
Biology 2022, 11, x FOR PEER REVIEW 3 of 17 to precipitate TFs, which are bound to their DNA targets. Therefore, the DNA sequences in the complex can be identified. An alternative technique has been applied with microarray technology [28], called protein-binding microarray (PBM), to identify TF binding sequences. Weirauch et al. [29] analyzed data from the PBM assays of over 1000 TFs to determine DNA binding specificity represented by position frequency matrices (PFMs). The PFMs were then used to infer PFMs of other TFs from diverse eukaryotes, including C. militaris, based on the similarity of protein sequences [29]. To construct a gene regulatory network (GRN), the Regulatory Sequence Analysis Tools (RSAT) allows users to scan the upstream sequences of genes with the PFMs using Markov chain-based background models to identify TF binding domains [30]. The predicted TF-target binding domain interactions from RSAT are then used to infer the GRN. This study aims to construct the GRN in C. militaris that responds to light. First, the collected PFMs of TFs of C. militaris were used to predict all possible TF binding sites in the promoter region of all genes in the C. militaris genome. Next, the interaction profiles between TFs and target genes were used to create a genome-scale GRN, which can be visualized as a network. Then, differentially expressed genes (DEGs) under light and dark conditions from transcriptome analysis were integrated into the network to extract a lightresponsive sub-network. The light-responsive GRN revealed the transcriptional regulation of many developmental and carotenoid and cordycepin biosynthesis genes by five key TFs, enabling more understanding of the regulation of light-controlled development and secondary metabolite production in C. militaris.

Materials and Methods
The overall workflow of the present study is shown in Figure 1 and is detailed as follows. Figure 1. The overall workflow implemented in the current study. First, position frequency matrices (PFMs) of each TF of C. militaris were collected from CIS-BP. Next, Matrix-scan, a package in the Regulatory Sequence Analysis Tools (RSAT), was used to scan the upstream sequences of all C. militaris genes with each PFM to identify TF binding sites for each TF. The TF-target gene interactions were used to construct a genome-scale GRN. Then, differentially expressed genes (DEGs) under light and dark conditions from transcriptome analysis were integrated into the network to extract a light-responsive GRN. Finally, a list of genes involved in the development and carotenoid and cordycepin biosynthesis were used to extract light-responsive regulatory sub-networks. The overall workflow implemented in the current study. First, position frequency matrices (PFMs) of each TF of C. militaris were collected from CIS-BP. Next, Matrix-scan, a package in the Regulatory Sequence Analysis Tools (RSAT), was used to scan the upstream sequences of all C. militaris genes with each PFM to identify TF binding sites for each TF. The TF-target gene interactions were used to construct a genome-scale GRN. Then, differentially expressed genes (DEGs) under light and dark conditions from transcriptome analysis were integrated into the network to extract a lightresponsive GRN. Finally, a list of genes involved in the development and carotenoid and cordycepin biosynthesis were used to extract light-responsive regulatory sub-networks.

Construction of the Genome-Scale GRN
The genome-scale GRN was constructed by identifying TF binding sites in the promoter region of all genes in the C. militaris genome. The genomic sequence of C. militaris was retrieved from BioProject (PRJNA41129) [31]. The upstream sequences from −1 to −1000 bp were extracted from 9651 genes. All 488 PFMs of 83 TFs of C. militaris were retrieved from the CIS-BP database version 2.0 (cisbp.ccbr.utoronto.ca, accessed on 6 December 2019) [29]. Matrix-scan, a package in the Regulatory Sequence Analysis Tools (RSAT) [30], was used to scan the upstream sequences with each PFM to identify TF binding sites for each TF in the genome. A higher-order (m = 2) Markov model was used as a background model with pseudo-frequencies = 0.001. Only predictions with a p-value ≤ 10 −5 were selected to reduce the false positives. The predicted interactions between TFs and their target genes were constructed as a genome-scale GRN in Cytoscape version 3.7.2 (https://cytoscape.org, accessed on 22 October 2019) [32], where nodes represent TFs and target genes and edges represent TF-target gene interactions.

Differentially Expressed Gene (DEG) Analysis and the Light-Responsive GRN
RNA sequencing reads of C. militaris strain TBRC6039 in light condition (BioProject: PRJNA579732, BioSample SAMN13118310) and dark condition (BioProject: PRJNA416937, BioSample SAMN07969453) were retrieved from previous studies [20,33]. Reads were filtered by the quality control at the Phred-score ≥30, yielding qualified reads of 88.68% and 94.46% in light and dark condition samples, respectively. The mapping process was performed using STAR two-pass-mode [34] version 2.7.3a. The reference genome of C. militaris was retrieved from BioProject (PRJNA41129) [31]. Uniquely mapped reads were more than 80%, while only 2.2% percent of the reads were unmapped in both light and dark conditions (Table S1 in Supplementary File S1). Non-uniquely mapped reads (18.4% and 16.1% of the total reads in light and dark conditions, respectively) and reads mapped with less than five reads (<0.005% of the total reads in both conditions) were removed. Finally, uniquely mapped reads with ≥5 mapped reads were counted by STAR using quantMode GeneCounts [34]. Read counts per gene in light and dark conditions were normalized using the trimmed mean of M values (TMM) method [35]. The normalization method is suitable for normalizing raw read counts between inter-and intra-samples before differentially expressed gene (DEG) analysis by NOISeq-sim [35]. The NOISeq-sim tool simulates the technical replication from the sequencing depth of input data [36]. Parameters of the simulation were set as follows: the percentage of the total number of reads (pnr) = 0.2, the number of simulated samples (nss) = 10 and variability in sample total reads of simulated samples (v) = 0.02. These parameters were set similarly to a previous analysis [20], except the value of nss, which was strictly set to reduce the false positives. NOISeq-sim calculates the probability of differential expression (q), which was set to be q ≥ 0.9 for DEGs in this study. The list of DEGs was then used to extract nodes and all connected edges among the nodes from the genome-scale GRN to construct the light-responsive GRN.

Identification of Key Regulator Genes
Key TFs that play an essential role in the global regulation of the network were identified based primarily on the outdegree of the nodes in the light-responsive GRN. Other node measures were also used for the analysis, including the normalized betweenness centrality score by cytoHubba, a plugin in Cytoscape [37] and the motif-based score calculated using an in-house script.
The betweenness centrality score of node n is defined as where σ ij is the number of shortest paths between node i and j and σ ij (n) is the number of shortest paths between node i and j that walk through node n. C represents the collection of all nodes in the network. The betweenness score of each node can be normalized by where B(n) is the betweenness score of node n, max(B) and min(B) are the maximum and minimum betweenness scores of all nodes in the network. The motif-based scores were calculated by counting the number of times a node participates in a feed-forward loop (FFL) motif as the master regulator, intermediate regulator, or target [38].

The Gene Regulatory Sub-Networks of Growth, Development, Carotenoid and Cordycepin Biosynthesis Pathways under Light Response
A list of genes involved in morphological development, carotenoid and cordycepin biosynthesis of model fungi was obtained from previous studies [12,[39][40][41][42][43][44] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/, accessed on 6 December 2019) [45]. The genes from C. militaris were identified based on homology analysis, which was done by bi-directional BLASTp with criteria evalue ≤10 −10 , percent identity ≥60% and percent coverage ≥40%. The development and biosynthesis subnetworks were created by extracting the light-responsive GRN with DEGs from the list and their parent TFs in the network.

The Genome-Scale GRN
The prediction results from RSAT's matrix-scan indicated 86,172 TF binding sites from 8584 genes. From a total of 83 TFs, 71 TFs were identified to regulate genes, resulting in 31,703 interactions between TFs and target genes ( Figure S1 in Supplementary File S2). The average outdegree of 71 TFs in the network was approximately 451, with a median of 161, indicating that some TFs regulate a large number of nodes in the network. All TF-target gene interactions in the genome-scale GRN are listed in Supplementary File S3.

The Light-Responsive GRN
To construct the light-responsive GRN, DEGs upon light stimulation were mapped into the genome-scale GRN and the DEG nodes along with edges among them were extracted. DEG nodes that did not connect to any other DEG nodes were removed. The extracted light-responsive GRN is visualized in Figure 2. Originally, the number of DEGs before mapping was 3256 (34.72% of the expressed genes), distributing into 1520 and 1736 up-regulated and down-regulated genes, respectively. After mapping and extracting, the light-responsive network contains 2689 genes (31.32% of nodes in the genome-scale GRN) and 6837 interactions (21.57% of interactions in the genome-scale GRN). There are 39 TFs out of 71 TFs (54.94%) in the network. TF-target gene interactions in the light-responsive GRN are listed in Supplementary File S4. Table 1 displays several node measures of 10 TFs in the light-responsive GRN. The 10 TFs are ranked in an order that achieves the maximum accumulated number of unique nodes directly regulated by the TFs (Figure 3). For example, the first TF, Snf21, an Snf2family helicase (CCM_04586), has the highest outdegree in the light-responsive GRN (1388 targets, approximately 52% of all nodes). The second TF is an AT-hook DNA-binding motif TF (CCM_08536) with an outdegree of 1044 target genes. Together, the AT-hook DNA-binding motif TF and Snf21 regulate 1861 unique target genes, approximately 70% of the total nodes. As shown in Figure 3, the top five TFs regulate about 90% (2396 genes) of the total nodes in the light-responsive GRN. Hereby, this study assigned these five TFs as the key regulators of the light-responsive GRN.
The five key TFs belong to general TF families identified in fungi [46]. In addition, the network reveals that the five key TFs regulate one another by forming several network motifs, as shown in Figure 4. Snf21 regulates a homeobox TF (CCM_07504), a forkhead box protein L2 (CCM_02646) and Heat shock factor 1 (Hsf1) (CCM_05142) in multiple feed-forward loops (FFLs). The AT-hook DNA-binding motif TF (CCM_08536) co-regulates the homeobox TF (CCM_07504) with Snf21. Finally, the homeobox TF and the forkhead box protein L2 engage in a feedback loop. Snf21 participates in 1407 feed-forward loops in the light-responsive network, where Snf21 is always the master regulator of the motifs (Table 1), implying a regulator role of Snf21 in the global control of the genes in C. militaris under light conditions. Further analysis shows that Snf21/Snf2 is an ortholog of SWI2 in baking yeast. Baking yeast's SWI2 is involved in the global transcriptional activation of many genes [47]. The protein also displays a conserved function for activating the transcription by controlling the chromatin structure [48,49]. In Schizosaccharomyces pombe, Snf2 is important to mitosis and vegetative growth [50]. TFs in the Snf2 family in Arabidopsis spp. and higher plants are involved in abiotic and oxidative stress responses [51][52][53].
Heat shock factor 1 (Hsf1) is another interesting TF in the light-responsive network. The TF has the highest betweenness score (a normalized betweenness score of 1.0), indicating it as the most central node that participates in the regulatory flow of other nodes in the network [38]. Hsf1 is activated by light and relates to fatty acid metabolism and circadian rhythm in eukaryotic organisms [54]. In yeast, Hsf1 is also required for growth at a normal temperature [55,56]. In C. militaris, besides light, heat is also a factor affecting cordycepin and carotenoid production [19]. The 10 TFs are ranked in an order that achieves the maximum accumulated number of unique nodes directly regulated by the TFs (see Figure 3). Outdegree represents the number of nodes directly regulated by the TF; Master means the number of times the TF acts as the master regulator of an FFL; Intermediate means the number of times the TF acts as the intermediate regulator of an FFL; Target means the number of times the node acts as the target of an FFL; Betweenness means normalized betweenness centrality score. Target means the number of times the node acts as the target of an FFL; Betweenness means normalized betweenness centrality score.

Growth and Developmental Regulation under Light Exposure
Growth and developmental genes in C. militaris were identified based on homology analysis against known genes in other model fungi (Table 2 and Supplementary File S5). To reveal how the genes are regulated, only DEGs in Table 2 and all upstream nodes from the light-responsive network were extracted. Figure 5 shows 15 growth and developmental genes that are included in the sub-network and are direct targets of the five key TFs.

Growth and Developmental Regulation under Light Exposure
Growth and developmental genes in C. militaris were identified based on homology analysis against known genes in other model fungi (Table 2 and Supplementary File S5). To reveal how the genes are regulated, only DEGs in Table 2 and all upstream nodes from the light-responsive network were extracted. Figure 5 shows 15 growth and developmental genes that are included in the sub-network and are direct targets of the five key TFs.   lreA (CCM_04514), lreB (CCM_00072), fphA (CCM_04461), veA (CCM_04531) and laeA (CCM_05395), encoding proteins that form the light-regulatory VeA/FphA/LreA/LreB complex [57] and the VeA/LaeA complex [42], were all differentially expressed upon light exposure. The complexes translocate into the nucleus and regulate asexual/sexual development balance [23]. VeA is a positive regulator of sexual fruiting body formation [58]. A veA mutation in A. nidulans reduced the development of sexual organs [23]. VeA has also been reported to be involved in the metabolism of nitrate and secondary metabolite in F. oxysporum [59].
Blue-light receptors LreA and LreB are activators of sexual development [57]. In contrast, red-light-mediated FphA inhibits fruiting body development while promoting conidiation in A. nidulans [42,60]. However, in C. militaris, light is essential for normal fruiting body formation [12]. LaeA is a negative regulator of sexual development [61] and VeA was shown to be the repressor of LaeA [59].
In the network ( Figure 5), Snf21 up-regulates lreB. Snf21 also forms an FFL with the homeobox TF (CCM_07504) to up-regulate and down-regulate veA and laeA, respectively. lreA was down-regulated in the transcriptomic data, but it was not included in Figure 5 because the analysis found no regulatory links between lreA and other genes. In addition, Hsf1 up-regulates the expression of fphA. Therefore, the analysis reveals Snf21, the homeobox TF (CCM_07504) and Hsf1 to be the main regulators of the VeA pathway in C. militaris. Evidence has supported that the homeobox TF (CCM_07504) was related to fruiting body formation and the mycelial growth of C. militaris [62]. Additionally, Snf21 up-regulates CCM_05639, importin β-2 subunit, which may involve in the translocation of the VeA complexes into the nucleus [63].
Snf21 and the homeobox TF (CCM_07504) also form an FFL to regulate oefC. The OefC protein promotes asexual development in A. nidulans [64]. The current study also identified fl (CCM_05556) as a target gene of OefC. In N. crassa, the fl gene encodes Fluffy, a transcriptional activator of several genes regulating conidiation [43]. In the network, both oefC and fl were down-regulated under light exposure. Consistent with a previous study, an oefC homolog in C. guangdongensis was also shown to down-regulate under light treatment [22].
In addition, Snf21 and Hsf1 down-regulate CCM_04849, which is identified as an flbC ortholog. Consistent with a previous study in C. militaris, it was found that CCM_04849 lreA (CCM_04514), lreB (CCM_00072), fphA (CCM_04461), veA (CCM_04531) and laeA (CCM_05395), encoding proteins that form the light-regulatory VeA/FphA/LreA/LreB complex [57] and the VeA/LaeA complex [42], were all differentially expressed upon light exposure. The complexes translocate into the nucleus and regulate asexual/sexual development balance [23]. VeA is a positive regulator of sexual fruiting body formation [58]. A veA mutation in A. nidulans reduced the development of sexual organs [23]. VeA has also been reported to be involved in the metabolism of nitrate and secondary metabolite in F. oxysporum [59].
Blue-light receptors LreA and LreB are activators of sexual development [57]. In contrast, red-light-mediated FphA inhibits fruiting body development while promoting conidiation in A. nidulans [42,60]. However, in C. militaris, light is essential for normal fruiting body formation [12]. LaeA is a negative regulator of sexual development [61] and VeA was shown to be the repressor of LaeA [59].
In the network ( Figure 5), Snf21 up-regulates lreB. Snf21 also forms an FFL with the homeobox TF (CCM_07504) to up-regulate and down-regulate veA and laeA, respectively. lreA was down-regulated in the transcriptomic data, but it was not included in Figure 5 because the analysis found no regulatory links between lreA and other genes. In addition, Hsf1 up-regulates the expression of fphA. Therefore, the analysis reveals Snf21, the homeobox TF (CCM_07504) and Hsf1 to be the main regulators of the VeA pathway in C. militaris. Evidence has supported that the homeobox TF (CCM_07504) was related to fruiting body formation and the mycelial growth of C. militaris [62]. Additionally, Snf21 up-regulates CCM_05639, importin β-2 subunit, which may involve in the translocation of the VeA complexes into the nucleus [63].
Snf21 and the homeobox TF (CCM_07504) also form an FFL to regulate oefC. The OefC protein promotes asexual development in A. nidulans [64]. The current study also identified fl (CCM_05556) as a target gene of OefC. In N. crassa, the fl gene encodes Fluffy, a transcriptional activator of several genes regulating conidiation [43]. In the network, both oefC and fl were down-regulated under light exposure. Consistent with a previous study, an oefC homolog in C. guangdongensis was also shown to down-regulate under light treatment [22].
In addition, Snf21 and Hsf1 down-regulate CCM_04849, which is identified as an flbC ortholog. Consistent with a previous study in C. militaris, it was found that CCM_04849 was regulated by light [12]. In A. nidulans, FlbC is required for proper conidiation initiation and the deletion of the gene delayed the conidiation process while enhancing fruiting body formation [65]. FlbC, along with Snf21 and the AT-hook DNA-binding motif TF (CCM_08536), regulates the expression of mat1 (CCM_06523), a sexual developmental gene [66]. FlbC and the AT-hook DNA-binding motif TF (CCM_08536) also regulate the expression of fadA. FadA is a signaling protein whose one of the substrate targets is PkaA, a regulator of sexual development [67]. In the network, the pkaA gene is up-regulated by Snf21.
Further, Hsf1 and the AT-hook DNA-binding motif TF (CCM_08536) down-regulate brlA (CCM_08959), an asexual developmental gene. Hsf1 also forms an FFL with the homeobox TF (CCM_07504) to up-regulate medA. BrlA was shown to transcriptionally regulate MedA and AbaA [44]. In the network, abaA is up-regulated by Snf21 and the AThook DNA-binding motif TF (CCM_08536). MedA and AbaA promote asexual growth [68] but also function in sexual development [44,69]. A medA mutant formed incomplete Hülle cells and did not produce cleistothecia [69]. Another asexual development-promoting gene identified in the network is flbB, which is down-regulated by the AT-hook DNA-binding motif TF (CCM_08536).
Taken together, the overall regulation in Figure 5 and expression profiles from the transcriptomic analysis suggest the transcriptional regulation of light that induces the growth and development of C. militaris into a sexual cycle by inactivating asexual developmental genes (flbB, flbC, oefC, brlA and fl) and activating sexual developmental genes (veA, fadA, pkaA and mat1).
It is important to note that the network in Figure 5 is far from complete. Many possible TFs were not assigned with the PFMs in CIS-BP and this study could not predict their regulatory links to target genes. For example, in A. nidulans, BrlA and AbaA are TFs that regulate one another [44,69]. These relationships were not identified in the present work.
The identified network in Figure 5 possesses several cascaded feed-forward loops (FFLs) made by regulatory interactions between the key TFs and other regulators. FFLs have been shown to play important roles in development, including temporal activation of genes or protection against transient loss of signals [70][71][72]. The current study identified Hsf1 as a target in two inter-connected FFLs, composed of Snf21-the homeobox TF-Hsf1 and Snf21-forkhead box L2-Hsf1. Based on the most frequent binding domain identified by RSAT (Figure 4), Snf21 binds to an AACAAT sequence, while the homeobox TF binds to TGACA, one of the conserved homeodomains reported in previous studies [73]. In addition, the forkhead box L2 recognizes the conserved domain GTAAACAA, similar to the forkhead box L1 [74]. Snf21 and Hsf1 form an FFL to regulate FlbC (CCM_04849). RSAT prediction identified (A)ACAAT and GAAnnTTCnnGAA domains in the promoter region of flbC. The (A)ACAAT is a TF-binding domain of Snf21, while GAAnnTTCnnGAA represents the heat shock element (HSE), a conserved binding domain of HSF with a repeated nGAAn or/and nTTCn sequence [75,76]. Furthermore, Snf21 and FlbC form another FFL to regulate mat1. Snf21 also forms several FFLs with the homeobox TF to regulate veA, laeA and oefC. We propose that these connected FFLs made of the key TFs and other regulators play important roles in the light-responsive regulation of growth and development in the fungi.

Light-Mediated Control of Secondary Metabolism
Light is an important factor affecting secondary metabolism in C. militaris. A study by Yang et al. [12] revealed that the fungi without the blue-light receptor gene (wc-1) showed a significant reduction in both cordycepin and carotenoid. Other evidence also suggested that the carotenoid content of C. militaris is enhanced by light [6,13], although light stress during the late growth stage may reduce carotenoid production. [19].
The carotenoid content measured in the culture from which the transcriptome in this study was extracted increased significantly in light compared to dark conditions [20]. Similar to a previous study [12], five orthologous genes of N. crassa carotenoid biosynthetic genes were not DEGs in this study (CCM_03697, CCM_03059, CCM_06355 geranylgeranyl diphosphate synthases, CCM_06728 carotenoid oxygenase 2 and CCM_09155 aldehyde dehydrogenase). Rather, it was identified that CCM_03203 FPP synthetase was up-regulated, while CCM_05579 geranylgeranyltransferase (GGTase) beta subunit was down-regulated. FPP synthetase catalyzes reactions that yield farnesyl pyrophosphate (FPP), a precursor of steroids and geranylgeranyl diphosphate (GGPP), a precursor of carotenoids [77]. GGTase is a CAAX protein prenyltransferase, utilizing GGPP as a substrate to perform the prenylation process on another CAAX-type protein, such as Ras1 [78,79]. Thus, the up-regulation of FPP synthetase and the down-regulation of GGTase may explain the increase of carotenoid production under light conditions [20]. In the light-responsive GRN, Hsf1 up-regulates the FPP synthetase gene (CCM_03203), while Snf21 down-regulates the GGTase gene (CCM_05579) (Figure 6). Similar to a previous study [12], five orthologous genes of N. crassa carotenoid biosynthetic genes were not DEGs in this study (CCM_03697, CCM_03059, CCM_06355 geranylgeranyl diphosphate synthases, CCM_06728 carotenoid oxygenase 2 and CCM_09155 aldehyde dehydrogenase). Rather, it was identified that CCM_03203 FPP synthetase was up-regulated, while CCM_05579 geranylgeranyltransferase (GGTase) beta subunit was down-regulated. FPP synthetase catalyzes reactions that yield farnesyl pyrophosphate (FPP), a precursor of steroids and geranylgeranyl diphosphate (GGPP), a precursor of carotenoids [77]. GGTase is a CAAX protein prenyltransferase, utilizing GGPP as a substrate to perform the prenylation process on another CAAX-type protein, such as Ras1 [78,79]. Thus, the up-regulation of FPP synthetase and the down-regulation of GGTase may explain the increase of carotenoid production under light conditions [20].
In the light-responsive GRN, Hsf1 up-regulates the FPP synthetase gene (CCM_03203), while Snf21 down-regulates the GGTase gene (CCM_05579) (Figure 6). Xia et al. [41] previously reported four enzymes (Cns1-Cns4) regulating the cordycepin production in C. militaris. Cns1/Cns2 is responsible for cordycepin biosynthesis. Cns3 and Cns4 are responsible for pentostatin (PTN) production and exportation, respectively. PTN is an adenosine deaminase inhibitor, which helps protect cordycepin from deamination. In the current study, cns1 (CCM_04436) was not a DEG, cns2 (CCM_04437) was down-regulated, while cns3 (CCM_04438) and cns4 (CCM_04439) were up-regulated. Thus, it is predicted that light may induce the production and exportation of PTN. As shown in Figure 6, cns2 is directly suppressed by Snf21 and the AT-hook DNA-binding motif TF (CCM_08536). cns3 is regulated by four TFs, including Snf21, the AT-hook DNA-binding motif TF (CCM_08536), Zn(II)2Cys6 TF OefC (CCM_05966) and FlbB (CCM_01128). OefC and FlbB are also regulators of asexual development [64]. Under light stimulation, both genes were down-regulated, allowing cns3 to be activated. oefC is suppressed by Snf21 and the homeobox TF (CCM_07504), while flbB is down-regulated by the AT-hook DNA-binding motif TF (CCM_08536) in the network ( Figure 6). OefC, along with the AT-hook DNA-binding motif TF and the homeobox TF, participates in the regulation of cns4. Taken together, we propose that PTN Xia et al. [41] previously reported four enzymes (Cns1-Cns4) regulating the cordycepin production in C. militaris. Cns1/Cns2 is responsible for cordycepin biosynthesis. Cns3 and Cns4 are responsible for pentostatin (PTN) production and exportation, respectively. PTN is an adenosine deaminase inhibitor, which helps protect cordycepin from deamination. In the current study, cns1 (CCM_04436) was not a DEG, cns2 (CCM_04437) was downregulated, while cns3 (CCM_04438) and cns4 (CCM_04439) were up-regulated. Thus, it is predicted that light may induce the production and exportation of PTN. As shown in Figure 6, cns2 is directly suppressed by Snf21 and the AT-hook DNA-binding motif TF (CCM_08536). cns3 is regulated by four TFs, including Snf21, the AT-hook DNA-binding motif TF (CCM_08536), Zn(II) 2 Cys 6 TF OefC (CCM_05966) and FlbB (CCM_01128). OefC and FlbB are also regulators of asexual development [64]. Under light stimulation, both genes were down-regulated, allowing cns3 to be activated. oefC is suppressed by Snf21 and the homeobox TF (CCM_07504), while flbB is down-regulated by the AT-hook DNAbinding motif TF (CCM_08536) in the network ( Figure 6). OefC, along with the AT-hook DNA-binding motif TF and the homeobox TF, participates in the regulation of cns4. Taken together, we propose that PTN production and an asexual cycle may be inversely coupled, where light and dark conditions promote the former and the latter, respectively.

Conclusions
This study combined the TF-target gene interactions predicted from RSAT with the transcriptomic analysis of C. militaris under a light-programming condition to construct the GRN that responds to light. In addition, the integration of the two information filters out false-positive predictions from both TF-target interaction and DEG analysis.
The current study identified five TFs (Snf21 (CCM_04586), an AT-hook DNA-binding motif TF (CCM_08536), a forkhead box protein L2 (CCM_02646), a homeobox TF (CCM_07504) and Hsf1 (CCM_05142)) that play key roles in the light-responsive regulation. The five TFs regulate almost over 90% of the total nodes in the light-responsive network, including genes involved in the growth and developmental programming and secondary metabolite production.
Additionally, this study identified a regulatory cascade made of several FFLs that co-regulate a large group of growth and developmental genes ( Figure 5). The analysis suggested that upon light stimulation, these connected FFLs regulate a decision switch between asexual and sexual growth by turning off several asexual developmental and conidiation-related genes, including flbB, flbC, oefC, brlA and fl while turning on genes involved in sexual development such as veA, fadA, pkaA and mat1. The light-mediated regulation is also coupled with genes involved in the biosynthesis of carotenoid and cordycepin/PTN ( Figure 6). Hsf1 up-regulates FPP synthetase (CCM_03203), while Snf21 down-regulates GGTase (CCM_05579), coupling fungal development with carotenoid production. Light also down-regulates OefC and FlbB, disabling asexual development while mediating the production of PTN, which may indirectly stabilize cordycepin from deamination. Figure 7 summarizes regulatory events that promote sexual growth and secondary metabolite production identified from the current study.
production and an asexual cycle may be inversely coupled, where light and dark conditions promote the former and the latter, respectively.

Conclusions
This study combined the TF-target gene interactions predicted from RSAT with the transcriptomic analysis of C. militaris under a light-programming condition to construct the GRN that responds to light. In addition, the integration of the two information filters out false-positive predictions from both TF-target interaction and DEG analysis.
The current study identified five TFs (Snf21 (CCM_04586), an AT-hook DNA-binding motif TF (CCM_08536), a forkhead box protein L2 (CCM_02646), a homeobox TF (CCM_07504) and Hsf1 (CCM_05142)) that play key roles in the light-responsive regulation. The five TFs regulate almost over 90% of the total nodes in the light-responsive network, including genes involved in the growth and developmental programming and secondary metabolite production.
Additionally, this study identified a regulatory cascade made of several FFLs that coregulate a large group of growth and developmental genes ( Figure 5). The analysis suggested that upon light stimulation, these connected FFLs regulate a decision switch between asexual and sexual growth by turning off several asexual developmental and conidiation-related genes, including flbB, flbC, oefC, brlA and fl while turning on genes involved in sexual development such as veA, fadA, pkaA and mat1. The light-mediated regulation is also coupled with genes involved in the biosynthesis of carotenoid and cordycepin/PTN ( Figure 6). Hsf1 up-regulates FPP synthetase (CCM_03203), while Snf21 down-regulates GGTase (CCM_05579), coupling fungal development with carotenoid production. Light also down-regulates OefC and FlbB, disabling asexual development while mediating the production of PTN, which may indirectly stabilize cordycepin from deamination. Figure 7 summarizes regulatory events that promote sexual growth and secondary metabolite production identified from the current study.  It should be noted that the other environment-influenced genes in our network still await further investigation. For example, several heat shock proteins were differentially expressed under light treatment, including CCM_06821, CCM_07508 and CCM_04804, which were also involved in the oxidative and osmotic stress responses in C. militaris in another study [80]. Additionally, two circadian clocks (CCM_01014; Frq and CCM_00908; Per) were down-regulated by the AT-hook DNA-binding motif TF (CCM_08536) and the homeobox TF (CCM_07504), respectively, under light exposure.
In summary, the construction of the light-responsive GRN reveals complex regulation between key TFs and target genes that leads to more understanding in light-regulating growth, development and secondary metabolite production in C. militaris. Nevertheless, further attempts to integrate the GRN with protein signaling and metabolic networks are needed to provide more complete descriptions of light-responsive processes in the fungi. This knowledge will enable more efficient approaches for the cultivation and beneficial metabolite production.