MYB4 is the best candidate transcription factor involved in pinosylvin stilbene biosynthesis in Pinus strobus L. cells by fungal elicitor treatment

Abstract

Key message

Dihydropinosylvin monomethyl ether (DPME) and pinosylvin monomethyl ether (PME) are pinosylvin derivatives that show high nematicidal activity against pine wood nematodes (PWNs). Here, we found that fungal elicitor treatment boosted the production of DPME and PME in cultured Pinus strobus L. cells and investigated the transcription factors (TFs) regulating the genes in the pinosylvin stilbenoid biosynthesis pathway.

Context

The discovery of TFs involved in the synthesis of DPME and PME provides an important clue to understanding the pinosylvin stilbenoid synthesis in pine plants.

Aims

We investigated the best fungal elicitor for the production of DPME and PME and the transcriptional activities of genes involved in PME and DPME biosynthesis in P. strobus L. cells after fungal elicitor treatment.

Methods

The content of DPME and PME in P. strobus cells was examined after treatment with fungal elicitors prepared from seven different species of fungi. Moreover, the role of fungal elicitors in the transcriptional activity of genes involved in DPME and PME biosynthesis was investigated by transcriptome analysis using RNA sequencing.

Results

Penicillium chrysogenum Thorn was the most efficient fungal elicitor for the production of DPME and PME in P. strobus cells among the other fungal species. The accumulation of DPME and PME in P. strobus cells after P. chrysogenum elicitor treatment increased 12.7-fold and 23.7-fold, respectively, compared to the control. Transcriptome analysis revealed that fungal elicitor treatment resulted in enhanced transcription of the PAL, 4CL, STS, PMT, and ACC genes, which are involved in PME and DPME biosynthesis. Some transcription factors belonging to the bHLH, MYB, WRKY, and ERF families showed a high transcription rate after fungal elicitor treatment.

Conclusions

We found that the fungal elicitor is a strong inducer of the accumulation of pinosylvin derivatives in the cells of P. strobus. We selected one unigene (c133966_g3_i1 in the MYB family) as the best candidate TF that regulates the DPME and PME biosynthesis in P. strobus by transcriptome analysis.

1 Introduction

Pine wilt disease, caused by the infection of pine wood nematode (PWN; Bursaphelenchus xylophilus) (Mamiya 1983; Naves et al. 2007; Futai 2013), seriously devastated the pine forests. Many Pinus species in East Asian countries and Europe are highly susceptible to PWN infection, but some species show resistance (Mamiya 1983; Hirata et al. 2017; Mariette et al. 2023). The detailed mechanism of PWN susceptibility and resistance in different pine species is not yet known. It has been suggested that defensive secondary compounds, preexisting and/or induced by PWN infection, may contribute to PWN resistance in some pine species (Suga et al. 1993; Hanawa et al. 2001; Hwang et al. 2021). The synthesis ability of secondary metabolites in response to PWNs can be related to the PWN resistance mechanism (Yamada and Ito 1993; Modesto et al. 2022).

Phytoalexins are important molecules in the context of plant immunity, as they play a crucial role in defending plants against various pathogens. Pinosylvin stilbenoids have antibacterial, antifungal, antifeedant, and nematicidal activities (Suga et al. 1993; Chong et al. 2009; Jeandet et al. 2010; Hwang et al. 2021). Pinosylvin and its methylated derivatives, pinosylvin monomethyl ether (PME) and dihydropinosylvin monomethyl ether (DPME), are representative phytoalexins found in pine species. Suga et al. (1993) reported that pinosylvin, DPME, and PME are the most toxic compounds against PWNs among the extracted compounds from the heartwood tissues of the Pinaceae family. Although pinosylvin and its methylated derivatives are very toxic to PWNs, these compounds are not detectable or exist in low amounts in fresh tissue of branches and needles of pine plants (Rosemann et al. 1991; Willför et al. 2003; Hwang et al. 2021) except for nonliving heartwood tissues (Harju and Venäläinen 2006). The production of pinosylvin stilbenoids in sapwood and needles in pine is only activated by abiotic and biotic stresses (Schoeppner and Kindl 1979; Gehlert et al. 1990; Rosemann et al. 1991; Harju et al. 2009).

Eastern white pine (Pinus strobus) displays resistance to PWNs (Mamiya 1983; Sutherland et al. 1991). Hwang et al. (2021) reported that PWN infection of P. strobus trees resulted in the strong enhancement of DPME and PME accumulation. They suggest that the ability of DPME and PME accumulation after infection of PWNs may be related to the PWN resistance of P. strobus. In contrast, PWN infection of two PWN-susceptible pine species (Pinus densiflora and Pinus koraiensis) did not show an increased accumulation of DPME and PME. The lack of PME and DPME accumulation after PWN infection in P. densiflora and P. koraiensis may be related to PWN susceptibility. This difference can be explained as acquired immune differences among the pine species against PWNs.

Plants have evolved immune systems to protect themselves from diseases or to defend themselves against pathogens such as bacteria, fungi, viruses, and pests. Understanding plant immunity, particularly for PWN resistance in pine plants, is essential for plant protection and development of PWN-resistant pine plants. Generally, there are two types of innate immunity against invading pathogens in plants: pathogen‐associated molecular pattern‐triggered immunity (PTI) and effector‐triggered immunity (ETI) (Tsuda and Somssich 2015; Yuan et al. 2021). PTI is triggered by the recognition of general pathogen patterns, while ETI is triggered by the recognition of specific pathogen effectors. Synthesis of phytoalexins is considered to be an important plant innate immune response to a variety of pathogens (Hain et al. 1993; Ahuja et al. 2012). Plant secondary compounds play an important role in preventing pathogen infection (Anjali et al. 2023). The phenylpropanoid pathway regulates plant basal immunity to viruses, bacteria, and fungi (Zabala et al. 2006; Konig et al. 2014; Zhang et al. 2020). The compounds of pinosylvin stilbenoids having a 1,2-diphenylethylene backbone are the groups in phenylpropanoid pathways.

Transcriptional regulation plays a key role in plant defense against pathogens (Li al. 2016). Transcription factors (TFs) play a crucial role in plant immunity by regulating the expression of genes involved in the plant’s defense responses against pathogens (Broun et al. 2006; Dubrovina and Kiselev 2017). TFs are essential components of both PTI and ETI in plant immunity (Tsuda and Somssich 2015). In both cases, TFs serve as central players in coordinating the plant’s defense mechanisms. They help regulate the expression of defense-related genes, allowing the plant to mount an effective response against invading pathogens. TFs are proteins that bind to specific DNA sequences in the promoter regions of target genes. In plant immunity, certain TFs are activated upon detection of pathogen-associated molecular patterns (PAMPs) or effectors produced by pathogens. Therefore, the discovery of TFs involved in the synthesis of pinosylvin derivatives provides an important clue to understanding the synthesis mechanism of pinosylvin derivatives in pine plants.

Several papers have identified TFs involved in the synthesis of the stilbene resveratrol, which is structurally similar to pinosylvin but has slightly different biosynthetic circuits. Resveratrol (3,4,5-trihydroxystilbene) is found in grapes, peanuts, and some berries of dicotyledonous plants (Tian and Liu 2020). Some of the TFs participating in regulating resveratrol biosynthesis have been reported in grapevine. R2R3-type MYB (v-myb avian myeloblastosis viral oncogene homolog) TFs, MYB14, and MYB15 were reported to regulate stilbene synthesis, increasing the promoter activity of Vitis STS genes (Höll et al. 2013). WRKY was also identified as the main TF family regulating resveratrol biosynthesis (Vannozzi et al. 2018). The conserved TF superfamily APETALA2/ethylene‐responsive factor (AP2/ERF) is involved in plant defense against pathogens, including bacteria, fungi, and nematodes (Moffat et al. 2012; Zheng et al. 2019; Zou et al. 2023). A Vitis quinquangularis TF, VqERF114, regulates stilbene synthesis by interacting with VqMYB35 (Wang and Wang 2019). The basic helix-loop-helix (bHLH) TF family is important for regulating the biosynthetic pathway of flavonoids such as flavones, flavonols, isoflavones, flavanones, flavanols, and anthocyanidins (Feller et al. 2011; Qian et al. 2021). There is little information on TFs regulating the synthesis of pinosylvin and its derivatives in pine species except for a report in which they suggested that a bHLH gene (PsbHLH1) strongly responds to PWN infection and may be related to the biosynthesis of pinosylvin stilbenoids in P. strobus (Hwang et al. 2021).

The enzymes pinosylvin synthase (STS) (Raiber et al. 1995) and pinosylvin O-methyltransferase (PMT) are critically involved in the biosynthesis of PME and DPME (Paasela et al. 2017). Transcription factors (TFs) play a key role in regulating their biosynthesis by controlling the expression of target genes (Broun et al. 2006; Dubrovina and Kiselev 2017). Therefore, the discovery of TFs involved in the synthesis of DPME and PME provides an important clue to understanding the synthesis mechanism of pinosylvin derivatives in pine plants.

Recently, Kim et al. (2022) reported that PME and DPME accumulation was highly increased in calli of P. koraiensis when the callus was treated with fungus-derived elicitors. In the cell culture of P. sylvestris, treatment with fungal mycelium extracts from Lophodermium seditiosum strongly increased the accumulation of PME hundreds of times (Lange et al. 1994). It is questioned whether fungal elicitor treatment increases PME and DPME accumulation in P. strobus cells and what kinds of TFs are involved in the accumulation of PME and DPME in P. strobus cells by fungal elicitor treatment.

Transcriptome analysis has been widely applied to examine the activity of genes in the biosynthesis of defense compounds in plants. To determine the response against PWN infection in pine plants, transcriptome analysis was carried out in P. massoniana (Liu et al. 2017), P. densiflora (Lee et al. 2019), and Pinus strobus plants (Hwang et al. 2021). There is no report on the transcriptional activity of genes involved in pinosylvin stilbenoid biosynthesis in pine plants after fungal elicitor treatment.

Here, we reported the time course accumulation of PME and DPME after fungal elicitor treatment in P. strobus cells. Comparative transcriptomic analysis of P. strobus cells was performed using the cDNA library from control and fungal elicitor-treated cells. We investigated the transcriptional activities of genes involved in pinosylvin-type stilbene biosynthesis and bHLH, MYB, WRKY, and ERF. A MYB gene was selected as the best candidate TF regulating the biosynthesis of DPME and PME by fungal elicitor treatment.

2 Materials and methods

2.1 Preparation of fungal elicitor

Penicillium chrysogenum (KCTC 6052), P. pinophilum (KCTC 16057), P. roquefortii (KCTC 6080), Trichoderma harzianum (KCTC 6426), Aspergillus niger (KCTC 6913), Aspergillus oryzae var. oryzae (KCTC 6983), and Aspergillus flavus (KCTC 6984) were purchased from the KCTC (Korean Collection for Type Cultures) in the Biological Resource Center of the Korea Research Institute of Bioscience and Biotechnology. The fungal cultures were maintained in potato dextrose agar (PDA) medium (Sigma, St. Louis, USA). The fungi were grown in 250-ml Erlenmeyer flasks containing 100 ml of potato dextrose broth (PDB) medium on gyratory shakers (110 rpm) for 10 days.

Fungal elicitors were prepared according to the procedure described by Baldi et al. (2009). The culture broths of fungi were obtained after 10 days of culture, and the culture broth was passed through a filter (Whatman No. 1). Then, the filtered medium was filtered again using a 0.22-µm filter for sterilization and used as a fungal elicitor.

2.2 Fungal elicitor treatment of P. strobus cells

Callus was obtained from mature seeds of P. strobus provided by the National Forest Seed and Variety Center in the Korea Forest Service, Chungju-si, Chungcheongbuk-do, 27,495, Republic of Korea. The protocol for callus induction followed the methods of a previous report (Koo et al. 2022). Calluses induced from cultured surfaces of zygotic embryos were subcultured at 3-week intervals. To establish the cell suspension culture, calli (500 mg) were transferred into 200-mL Erlenmeyer flasks containing 50 mL 1/2 LV liquid medium (Litvay et al. 1985) with 20 g/L sucrose with 0.5 mg/L 2,4-D and 0.5 mg/L BA (benzyl adenin). The suspension culture was maintained by subculture at 10-day intervals. The flasks were shaken at 120 rpm in the dark.

To investigate the best fungal elicitor among the prepared seven different fungi mentioned above, the fungal elicitors were treated in culture medium at a concentration of 2.0% (v/v) on 10 days of cell suspension culture. The accumulation of DPME and PME was examined after 3 days of fungal elicitor treatment. The control received an equal amount of PDA liquid medium instead of a fungal elicitor solution. Cells were harvested after 3 days of culture, and the contents of PME and DPME were analyzed by GC/MS.

To investigate the time-dependent accumulation of DPME and PME in cultured cells of P. strobus after fungal elicitor (P. chrysogenus) treatment, the cells were harvested after 0, 1, 3, and 5 days of culture, and the contents of PME and DPME were analyzed by GC/MS. All experiments were performed three times with triplicate samples.

2.3 Analysis of DPME and PME in cells by GC‒MS

Sampled cells for GC analysis were dried at 50 °C in an oven. The powder (200 mg), prepared by grinding each sample, was soaked in 100% methanol (1 ml) and subjected to sonication for 30 min at 40 °C. After centrifugation (15,000 × g for 10 min), the supernatant was filtered using a syringe filter. The extract aliquots were analyzed by gas chromatography (GC, Agilent 7890A) and mass spectrometry (MS, Agilent 5975C) equipped with an HP-5MS capillary column. The injection temperature was 250 °C with a split/splitless injection (10:1), and the column temperature program was 70 °C for 4 min, followed by an increase to 220 °C at 5 °C/min. The interface temperature of GC/MS was set at 300 °C. The temperature of the ionization source chamber was 250 °C. The standards of PME and DPME compounds were purchased from Sigma–Aldrich Co.

2.4 Preparation of cDNA library and RNA sequencing

Total RNA was extracted from the cultured cells treated with and without fungal elicitor using the RNeasy® Plant Mini Kit (Qiagen, Germany). DNA contamination was eliminated from the total RNA using DNase. mRNA was isolated from total RNA using an mRNA purification kit (Stratagene, USA). Purified mRNA was fragmented by the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, CA, USA) for short-read sequencing. The fragmented RNAs were reverse-transcribed into cDNA. Then, ligate adapters were applied to both ends of the cDNA fragments. After amplifying fragments using PCR, fragments with insert sizes between 200 and 400 bp were selected. The paired ends of the libraries were sequenced with a 100-bp read length using Illumina sequencing by MacroGen Inc. Seoul, Korea (www.macrogen.com).

2.5 De novo transcriptome assembly and differentially expressed gene analysis

With these low sequence data, sequence trimming was performed using the sliding window method with the Trimmomatic program (Bolger et al. 2014) (Appendix Table 2). The trimmed reads are merged into one file to perform transcriptome assembly. Merged data were assembled using the Trinity program (Grabherr et al. 2011), and the combined read sequences were converted into contigs. Read mapping was performed to identify unigenes using Bowtie software (Langmead et al. 2009). The assembled contigs were filtered and clustered using the CD-HIT program to obtain nonredundant transcripts (Li and Godzik 2006). Assembled contigs are represented as expressed transcripts. Contigs were merged into nonredundant unique transcripts and clustered into unigenes.

The abundance of unigenes across samples is estimated by the RSEM algorithm (Li and Dewey 2011). The expression value calculated by read count was used for further differentially expressed gene (DEG) analysis. DEG analysis was performed on a comparison pair (Fungi vs Wt) as requested using edgeR. The results showed 12,910 contigs that satisfied arbitrary fold change (FC) values > 2 and a significance level set at p values of < 0.05. Multiplicity correction was performed by applying the Benjamini‒Hochberg method (1995) to the p values to control the false discovery rate (FDR).

2.6 Functional annotation and enrichment analysis

To annotate the clustered unigenes, they were blasted against public databases with BLASTN and BLASTX, including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Pfam, NCBI Nucleotide (NT), NCBI nonredundant Protein (NR), UniProt and EggNOG, using DIAMOND software with an E-value default cut-off of 1.0E − 5 (Buchfink et al. 2021). The unigenes were also mapped back to the Cluster of Orthologous Groups of Proteins (COG) database to analyze functional categories.

To conduct functional enrichment analysis, we counted the number of total genes and DEGs for each KEGG pathway using homebrew Python code. Then, the statistical significance was determined by hypergeometric tests using the Python scipy module. The 30 most significantly overrepresented pathways were visualized using the Python seaborn module.

2.7 qPCR analysis of selected genes

To determine the effect of fungal elicitor treatment on the transcriptional activities of stilbene biosynthesis and their related TF genes, P. strobus cells were treated with or without fungal elicitor, and the control was treated with distilled water. The cultured cells were collected at 0, 3, and 7 days after fungal elicitor treatment and immediately stored at − 80 °C until use.

RNA was extracted from control and fungal elicitor-treated cells and reverse transcribed using M-MLV Reverse Transcriptase (Invitrogen, CA, USA). The Rotor-Gene SYBR Green PCR Kit (Qiagen, Germany) was used for qPCR analysis. The qPCR analysis was performed with three biological replicates and presented as the average relative quantities ± SE. The relative expression value for each gene was calculated using the 2^−ΔΔCT method (Livak and Schmittgen 2001). The β-actin gene of P. strobus was used for normalization. The primers used in this study are listed in Appendix Table 3.

2.8 Protein network analysis

Protein network analysis was performed to evaluate the degree of connectivity between the various TFs and proteins that participate in pinosylvin biosynthesis. The analysis was performed using STRING v11.0 (https://string-db.org/) with a minimum needed interaction score set at “medium confidence” (0.3). Nine TF proteins and five proteins involved in pinosylvin biosynthesis in P. stroubus were chosen as query proteins to search the STRING database. To identify TF regulators of DEGs involved in PME and DPME biosynthetic pathways, an integrative regulatory network analysis was conducted by iRegNet protocol (Shim et al. 2021) for Arabidopsis homologs of five DEGs in the stilbenoids, diarylheptanoids, and gingerol biosynthetic pathways (ko00945) and five DEGs coding PAL, 4CL, STS, and PMT, respectively.

2.9 Statistical analysis

Analysis of samples for GC/MS and qPCR analysis were performed in triplicates. Values in all data are presented as the average relative quantities ± standard error (SE). Statistical analysis was performed using SPSS software (SPSS Science, Chicago, IL, USA).

3 Results

3.1 DPME and PME accumulation in P. strobus cells after fungal elicitor treatment

Seven different fungal elicitors were added to P. strobus cells at 10 days of culture. The concentration of the fungal elicitor was adjusted to 2% v/v in the liquid cell culture medium. The treatment of all fungal elicitors prepared from seven species of fungi resulted in a high increase in DPME and PME accumulation in cells after 3 days compared to the control (Fig. 1a–h). GC total ion chromatogram revealed that fungal elicitor treatment particularly stimulated the accumulation of DPME and PME. Their accumulation differed greatly depending on the fungal species (Fig. 1j–k). The DPME and PME in P. strobus cells were identified by comparison with standard compounds (Fig. 1i). Overall, Penicillium-derived elicitors were superior to other species, and in particular, the elicitor made from P. chrysogenum was the best. Trichoderma and Aspergillus species showed lower elicitor effects than Penicillium species. Therefore, future studies should use a fungal elicitor derived from P. chrysogenum.

Fig. 1

GC/MS analysis of DPME and PME accumulation in cultured cells of P. strobus by seven different fungal elicitor treatments. a GC total ion chromatogram of methanol extracts of control cells without elicitor treatment. b–h Chromatograms of methanol extracts of cells after seven different fungal elicitor treatments. i Chromatogram of authentic standards of DPME and PME. j Content of DPME in cells after 3 days of fungal elicitor treatments. k Content of PMEs in cells after 3 days of fungal elicitor treatments

The accumulation of DPME and PME in P. strobus cells was monitored at 0, 1, 3, and 5 days of P. chrysogenum elicitor treatment. The accumulation of DPME and PME was significantly increased by elicitor treatment and showed a peak at 3 days (Fig. 2). In the control cells without elicitor treatment, the contents of DPME and PME in cells were 161 ± 12 µg/g DW and 11.2 ± 1.7 µg/g DW, respectively (Fig. 2). At 3 days of elicitor treatment, the DPME content showed a 12.7-fold increase (2047 µg/g DW), and the PME content showed a 23.7-fold increase (265 µg/g DW) compared to the control without elicitor treatment (Fig. 2). In all treatments, the content of DPME was higher than that of PME (Fig. 2).

Fig. 2

Content of DPME and PME in cells during 5 days of culture after elicitor treatment. Different letters above the bars indicate significantly different values (P < 0.05), calculated using one-way ANOVA followed by Duncan’s post hoc analysis

3.2 De novoassembly of the P. strobus transcriptome

It was speculated that the enhanced production of DPME and PME in the cultured cells of P. strobus by fungal elicitor treatment is associated with altered gene expression in the pinosylvin stilbenoid biosynthesis pathway. To elucidate the molecular mechanism of the enhanced production of DPME and PME biosynthesis by fungal elicitor treatment, mRNA was extracted from cultured cells with or without elicitor treatment, and the transcripts were sequenced using an Illumina platform.

Read mapping was performed to identify unigenes using the trimmed read information for each sample using Bowtie software. As a result, the mapping rate for each sample was 70.34 ~ 79.75% (Appendix Table 2). The total number of transcripts obtained from the six samples was 201,946, and the total number of genes was 155,625, with the average length of each sequence being 363 bp. With a total of 155,625 transcripts, 135,735 unigenes were obtained after filtering the longest contig. These unigenes were used for analysis.

As a result of performing alignment with reference genes using various analysis programs for functional annotation analysis, 45,302 out of a total of 135,735 unigenes were annotated, showing an annotation rate of 33.38% (Appendix Fig. 12). Boxplots show the expression distribution of read counts between control and fungal elicitor-treated samples based on raw data, log2 transformation of reads, and TMM normalization (Appendix Fig. 13a). The density plot shows the expression distribution of the entire transcriptome for each sample to raw data, log2 transformation of reads, and TMM normalization (Appendix Fig. 13b).

3.3 Identification of DEGs

Among a total of 135,735 contigs obtained from 6 samples, low-quality transcripts were filtered and then TMM normalized to select 31,971 contigs. Furthermore, a total of 12,910 DEGs belonging to log2 (fold change (FC) ≥ 2 and p value < 0.05) were finally selected. Through DEG analysis, 6589 DEGs were downregulated and 6321 DEGs were upregulated in the fungal elicitor treatment group compared to the control group (Appendix Fig. 13c). To visualize significant differences in DEGs, a volcano plot (Appendix Fig. 13d) and smear plot (Appendix Fig. 13e) were prepared. Hierarchical clustering heatmap analysis showed a large difference in expression between the control group and the elicitor treatment group (Fig. 3a), which is consistent with the multidimensional scaling analysis (Fig. 3b) that grouped control and fungal elicitor treatment in different axes.

Fig. 3

Hierarchical clustering and multidimensional scaling analysis results. a Heatmap result of hierarchical clustering analysis, which clusters the similarity of contigs and samples by expression level (normalized value). b Multidimensional scaling analysis between control and fungal elicitor-treated samples

3.4 Functional classification of DEGs

As a result of Gene Ontology (GO) analysis and COG functional classification with DEG for functional annotation of Unigenes, biological processes (BP), cellular component (CC), and molecular function (MF) results are shown (Fig. 4). As a result of BP analysis, the expression of genes related to cellular processes and metabolic processes related to the production of secondary metabolites was high in cultured cells (Fig. 4a). The expression of genes was the highest in the cell part in CC and catalytic activity in the MF category (Fig. 4a). In COG functional classification, the cluster of “replication, recombination and repair” showed the largest proportion of the group, followed by “signal transduction mechanism” and “posttranslational modification, protein turnover, chaperones” (Fig. 4b).

Fig. 4

GO analysis and COG classification results of the P. strobus cell transcriptome. A total of 12,910 unigenes were annotated using DIAMOND software and summarized into the three GO categories. a Diagrams of GO analysis. b COG functional classification

To gain deeper insight into the biological functions of DEGs, we performed functional enrichment analysis on 12,910 DEGs based on KEGG annotations. Interestingly, molecular pathways related to secondary metabolite biosynthesis (ko01110 and ko00999) and phenylpropanoid biosynthesis (ko00940) were overrepresented (Fig. 5). Furthermore, three out of five genes involved in the stilbenoid, diarylheptanoid, and gingerol biosynthetic pathways were transcriptionally controlled (enrichment analysis p value < 0.01). These results suggested that genes involved in pinosylvin biosynthesis might be affected by treatment with fungal elicitors. Thus, we further investigated the transcriptional activity of genes involved in PME and DPME biosynthesis.

Fig. 5

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for upregulated DEGs in P. strobus cell transcriptomes. Functional enrichment analysis was conducted on the same set of 12,910 unigenes based on the KEGG annotations. The x and y axes represent the overexpressed pathways. The color of the dots indicates a higher-order classification of the enriched pathway. The size of the dots represents the number of DEGs to the total number of DEGs in a particular pathway

3.5 Transcriptional activities of genes involved in the PME and DPME biosynthetic pathways

PME and DPME biosynthesis started from the enzymes phenylalanine ammonia-lyase (PAL), followed by 4-coumarate-CoA ligase (4CL), acetyl-CoA carboxylase (ACC), pinosylvin synthase (STS), and pinosylvin O-methyltransferase (PMT), as shown in Fig. 5. Contigs for the putative PAL, 4CL, ACC, STS, and PMT genes were selected. It was confirmed that PAL, 4CL, STS, and PMT contigs exist as a family except for the ACC contig. The FPKM values were different among the contigs. We selected each contig showing highly enhanced transcription values after fungal elicitor treatment, and c46784_gi_i1 for PAL, c132355_g8_i3 for 4CL, c143977_gi_i1 for STS, and c134134_gi_i1 for PMT were selected for further qPCR analysis (Fig. 6). Pinosylvin is produced by the condensation of cinnamoyl-CoA with three molecules of malonyl-CoA. ACC is involved in the conversion of acetyl-CoA to malonyl-CoA (Sasaki and Nagano 2004). One homomeric ACC contig (c116421_g1_i1) was found in the annotated transcriptomes, and the expression of the gene was increased by fungal elicitor treatment (Fig. 6).

Fig. 6

Biosynthetic pathway of DPME and PME and heatmap analysis for the expression of differentially expressed genes (DEGs) associated with DPEM and PME biosynthesis between control (C1–C3) and fungal elicitor-treated samples (T1–T3). PAL phenylalanine ammonia lyase, 4CL 4-coumaroyl-CoA ligase, STS stilbene synthase, PMT pinosylvin O-methyltransferase, ACC acetyl-CoA carboxylase

3.6 Identification of transcription factors

bHLH, MYB, WRKY, and ERF contigs were selected from the transcriptome sequences obtained from cultured P. strobus cells. As a result of heatmap analysis, a total of six contigs among bHLH, MYB, WRKY, and ERF contigs showed the highest increase in FPKM value in the fungal elicitor treatment compared to the control group (Fig. 7, arrows). The two bHLH TFs (c126595_g1_i1 and c161886_g1_i1) had helix–loop–helix motifs and were related to functionally unknown TF genes of Arabidopsis bHLH144 and Arabidopsis bHLH147, respectively (Table 1). The two WRKY TFs (c48708_g1_i1 and c63840_g1_i1) are similar to Arabidopsis WRKY75 and WRKY6, both of which are defense-related TFs (Table 1). The two c52919_g1_i1 and c133966_g3_i1 TFs had R2R3-MYB motifs and were related to Arabidopsis MYB112 and Arabidopsis MYB4, both of which are involved in flavonoid biosynthesis (Table 1). The two ERF TFs, c130895_g1_i1 and c141277_g1_i1, were related to Arabidopsis ORA59, Arabidopsis ERF2, and Arabidopsis ERF110; all of which are involved in the ethylene signaling network (Table 1).

Fig. 7

Heatmap analysis for the expression of differentially expressed genes (DEGs) in bHLH (a), WRKY (b), MYB (c), and ERF (d) TFs between control (C1–C3) and fungal elicitor-treated samples (T1–T3). Arrows indicate the contigs showing highly enhanced FPKM values after fungal elicitor treatment compared to the control

Table 1 Selection of contigs showing highly enhanced read number in bHLH, WRKY, MYB, and ERF TFs after fungal elicitor treatment in annotated unigenes of the P. strobus transcriptome

3.7 qPCR analysis of stilbene biosynthetic pathway genes and selected TF genes

qPCR analysis revealed that the expression levels of selected PAL, 4CL, ACC, STS, and PMT contigs were enhanced after fungal elicitor treatment (Fig. 8a–d). STS showed the highest expression after fungal elicitor treatment among the four tested genes. PAL, 4CL, ACC, and STS showed the highest expression at 7 days of fungal elicitor treatment. PMT showed the highest expression at 3 days of fungal elicitor treatment and declined at 7 days (Fig. 8e).

Fig. 8

Transcriptional activity of stilbene biosynthetic pathway genes and selected TF genes. Expression of the PAL (a), 4CL (b), STS (c), PMT (d), and ACC (e) genes in P. strobus cells after fungal elicitor treatment by qPCR analysis. f Transcriptional activity of two bHLH contigs (c126595_g1_i1 and c161886_g1_i1). g Transcriptional activity of two WRKY contigs (c48707_g1_i1 and c63840_g1_i1). h Transcriptional activity of two MYB contigs (c52919_g1_i1 and c133966_g3_i1). i Transcriptional activity of two MYB contigs (c130895, c134931, and c141277). The expression data were normalized to β-actin. The vertical bars in each column represent the mean ± SE based on three biological replicates

A total of nine sequences of TFs showed highly increased FPKM values after fungal elicitor treatment by qPCR analysis (Fig. 8f-i). The efficiency of qPCR reactions for selected genes is essential for obtaining accurate and reliable quantitative data. qPCR efficiency of nine selected TFs revealed that the qPCR efficiency of all the 10 genes lies between 90 and 110% efficiency (Appendix Fig. 14). All six unigenes of bHLH, R2R3-MYB, and WRKY TFs were highly enhanced after fungal elicitor treatment (Fig. 8f-h). Two MYB unigenes (c52919_g1_i1 and c133966_g3_i1) were more quickly enhanced at 3 days of fungal elicitor treatment compared to other bHLH and WRKY contigs. The unigenes of bHLH and WRKY TFs showed the highest expression at 7 days of fungal elicitor treatment. Among the three unigenes of ERF TFs, c130895_g1_i1 showed the highest expression after fungal elicitor treatment compared to the other two unigenes (c134931_g1_i1 and c141277_g1_i1) (Fig. 8i).

3.8 Identification of TFs by network analysis with DEGs involved in PME and DPME biosynthetic pathways

To identify potential master regulators of DEGs involved in PME and DPME biosynthetic pathways, we conducted by iRegNet protocol (Shim et al. 2021) for Arabidopsis homologs of five DEGs in the stilbenoid, diarylheptanoid, and gingerol biosynthetic pathways (ko00945) and five DEGs coding PAL, 4CL, STS, and PMT, respectively. Notably, our analysis revealed multiple MYB family transcription factors (TFs) as putative upstream regulators of these genes (Fig. 9a–b). Additionally, our investigation revealed a close association between MYB4 and multiple PAL and 4CL in Arabidopsis (Fig. 9c).

Fig. 9

Putative upstream regulatory transcription factors of DEGs involved in PME and DPME biosynthetic pathways depicted using iRegNet analysis (Shim et al. 2021). a High-confidence upstream transcription factors of regulating Arabidopsis homologs of five DEGs in the stilbenoid, diarylheptanoid, and gingerol biosynthetic pathways. b High-confidence upstream transcription factors of regulating Arabidopsis homologs of five DEGs coding PAL, 4CL, STS, and PMT. Orange and Blue circles indicate DEGs and upstream TFs, respectively. c Molecular connections between MYB4, PAL, and 4CL. Molecular links of MYB4 were retrieved from STRING database

3.9 Network analysis of genes of the stilbene biosynthetic pathway and selected TFs

We performed a correlation analysis between the genes involved in pinosylvin biosynthesis and significant TFs identified in fungal elicitor-treated groups based on the protein interaction network obtained from STRING (Fig. 10). Out of 14 submitted proteins, the database identified 11 proteins, and the interaction analyses were shown for these 11 protein nodes having 18 edges with 3.27 average node degree, 0.73 average local clustering coefficient, and a PPI enrichment p value of 6.54e − 13. MYB4 (P. strobus contig, c133966_g3_i1) was connected with 4CL2, OMT1, CHS, and PAL2. Moreover, MYB4 was also connected to TFs (MYB78 and WRKY75). WRKY75 showed a central interconnect among the TFs because the gene was connected to MYB4, MYB78, ERF110, and WRKY6. However, two bHLHs (P. strobus contig, c126595_g1_i1 and c161886_g1_i1) and ERF (P. strobus contig, c141277_g1_i1) showed no homology to the protein of Arabidopsis.

Fig. 10

Predicted protein‒protein interaction network between the genes involved in pinosylvin biosynthesis and significant TFs identified in fungal elicitor-treated groups. Arabidopsis thaliana was used as the reference to obtain homologous proteins. Edges represent protein‒protein associations. The analysis was performed using STRING v11.0 with the minimum needed interaction score set at “medium confidence” (0.3)

4 Discussion

Pinosylvin stilbenoids are known as representative phytoalexins in pine species and have strong nematicidal activities against PWNs (Suga et al. 1993; Hwang et al. 2021). Thus, it is expected that pinosylvin stilbenoids can be developed into eco-friendly PWN prevention materials. In the current study, fungal elicitor treatment of P. strobus cells also resulted in a high accumulation of DPME and PME. These results indicate that the fungal elicitor treatment may be directly related to the stimulation of biosynthesis and accumulation of pinosylvin derivatives in P. strobus cells.

Among the fungal elicitors derived from seven different fungi, a fungal elicitor prepared from P. chrysogenum was the most effective elicitor for boosting the production of DPME and PME in P. strobus cells. The genus Penicillium comprises ascomycetous fungal species with important implications for natural environments as well as in the food and drug industries. In particular, P. chrysogenum is the main source of penicillin (Barreiro and García-Estrada 2019).

Recently, there have been great advances in fungal elicitor-triggered plant immunity, which provides insights into plant‒microbe interactions and plant disease control (Guo and Cheng 2022). The molecular defense response of pine trees (Pinus spp.) by elicitor treatment to control PWNs was reviewed by Modesto et al. (2022). The application of chitosan elicitor to pine trees before PWN infection has been investigated as a PWD control method (Nunes da Silva et al. 2021). The addition of biofertilizer (mixed with diazotrophic bacteria and a chitosan-producing fungus, Cunninghamella elegans) to soils confers resistance against the PWN in P. pinaster (Nunes da Silva et al. 2019). Elicitation treatment using diketopiperazine derived from Bacillus thuringiensis can control pine wilt disease (Park et al. 2020). These results suggest that elicitors can be used as economical and eco-friendly biocontrol agents to prevent pine wilt disease in PWN-susceptible pines.

We performed transcriptome analysis to understand the expression mechanisms of genes involved in the enhanced production of DPME and PME in cultured cells of P. strobus by elicitor treatment. The expression levels of genes involved in the synthesis of pinosylvin stilbenoids were compared by dividing the control group and the fungal elicitor treatment group. The PAL, 4CL, STS, PMT, and ACC genes are involved in the biosynthesis of PME and DPME from phenylalanine (Dubrovina and Kiselev 2017). Among the annotated unigenes, the expression level was highly increased in 1 (c46784_g1_i1) out of six PAL contigs and one (c132355_g8_i3) out of eight 4CL contigs, one (c143977_g1_i1) out of two STS contigs, and one (c134134_g1_i1) out of seven PMT contigs. Previously, we observed the expression patterns of PAL, 4CL, STS, and PMT unigenes in P. strobus plants after PWN infection (Hwang et al. 2021). These four contigs obtained in transcriptome analysis of P. strobus cell were exactly matched with the selected four contigs in the transcriptome data of P. strobus after PWN infection. The four selected PAL, 4CL, STS, and PMT unigenes may have primarily participated in the biosynthesis of pinosylvin derivatives regardless of PWN infection or fungal elicitor treatment.

Transcription factors are essential components of both PTI and ETI in plant immunity (Tsuda and Somssich 2015). Transcriptional reprogramming regulated by TFs occurs by immediate responses to pathogen infection. In general, TFs bind to the promoters of genes to regulate the transcription of genes (Roeder 1996; Latchman 1997). In the biosynthesis of resveratrol stilbene, particularly in grapes, few TFs related to the synthesis of resveratrol stilbene have been reported in grapes. R2R3-type MYB TFs (VviMYB14 and VviMYB15) have been reported as representative transcriptional regulators involved in resveratrol synthesis in grapes, and these two genes affect the promoter activity of the STS gene (Höll et al. 2013). Moreover, WRKY family transcriptional regulators also regulate resveratrol synthesis in grapes (Wang et al. 2020; Vannozzi et al. 2018). There has been no report on the direct role of the AP2/ERF gene in the biosynthesis of stilbenes. However, various AP2/ERF TF-binding sites are present within the promoter region of the grapevine STS gene (Wang and Wang 2019).

Until now, TFs regulating the synthesis of pinosylvin and its derivatives in pine species have not yet been identified. In the present study, we examined the expression profile of bHLH, R2R3-type MYB, and WRKY and ERF TFs in P. strobus cells after fungal elicitor treatment. Because the transcriptome data of P. strobus plants after infection of PWNs are available (Hwang et al. 2021), we tried to compare the expression patterns of the bHLH, R2R3-type MYB, WRKY, and ERF TF unigenes in the transcriptome data between PWN infection and fungal elicitor treatment. We selected several unigenes in the family of bHLHs (c126595_g1_i1 and c161886_g1_i1), R2R3-type MYBs (c52919_g1_i1 and c133966_g3_i1), and WRKY TFs (c48708_g1_i1 and c63840_g1_i1), which had highly increased FPKM values in the transcriptome of P. strobus cells after fungal elicitor treatment. However, the selected unigenes in the transcriptome data by fungal elicitor treatment did not show the same pattern in the transcriptome of P. strobus plants after PWN infection except for one unigene (c48708_g1_i1) of the WRKY TF (Hwang et al. 2021). Interestingly, the unigene (c48708_g1_i1) is highly expressed in P. strobus not only after PWN infection but also after fungal elicitor treatment. The deduced amino acid sequences of c48708_g1_i1 had 78% identity to those of Arabidopsis WRKY75, which positively regulates jasmonate (JA)-mediated plant defense against necrotrophic fungal pathogens (Chen et al. 2021).

In the finding of master regulators of DEGs involved in PME and DPME biosynthetic pathways, iRegNet analysis (Shim et al. 2021) was conducted using Arabidopsis homologs of five DEGs in the stilbenoid, diarylheptanoid, and gingerol biosynthetic pathways (ko00945) and five DEGs coding PAL, 4CL, STS, and PMT, respectively. This analysis revealed that MYB family TFs were the main putative upstream regulators of these genes, suggesting a potentially pivotal role for MYB family TFs in stilbenoid biosynthesis. Focused analysis on MYB TFs within the DEGs revealed that two up-regulated MYB TFs, namely MYB112 (c52919_g1_i1) and MYB4 (c133966_g3_i1), are strongly responsive following fungal elicitor treatment. The expression patterns of two MYB unigenes (c52919_g1_i1 and c133966_g3_i1) showed rapid enhancement at 3 days of fungal elicitor treatment. Additionally, MYB4 showed a close association among multiple PAL and 4CL in Arabidopsis. By network analysis of genes by String, the contig of c133966_g3_i1 homologous to Arabidopsis MYB4 was well connected with 4CL2, OMT1, CHS, and PAL2. Moreover, MYB4 was also connected to other TFs (MYB78 and WRKY75). Although WRKY75 was not directly connected to the genes in pinosylvin stilbenoid biosynthesis, WRKY75 showed the highest connection among TFs. Our finding strongly suggests the potential role of MYB4 (c133966_g3_i1) as a master regulator in the biosynthesis of PME and DPME in response to fungal elicitor treatment.

Figure 11 shows the simplified representation of the production of nematicidal compounds (PME and DPME) in the culture of P. strobus cells by fungal elicitor treatment. A candidate TF gene putatively involved in PME and DPME biosynthesis was selected by transcriptome analysis of P. strobus cells.

Fig. 11

Simplified representation for both the analysis of nematicidal compounds (PME and DPME) in P. strobus cells by fungal elicitor treatment and identification of TF genes involved PME and DPME biosynthesis by transcriptome analysis

5 Conclusion

In this experiment, treatment of P. strobus cells with fungal elicitors significantly increased the synthesis of pinosylvin stilbenes. This result can be interpreted as a result of the increased expression of genes involved in the synthesis of pinosylvin stilbenes by the elicitor treatment. In general, the synthesis of secondary metabolites in plants is known to be centrally regulated by TFs. Therefore, to select TFs that control pinosylvin stilbene synthesis upon fungal elicitor treatment, we treated P. strobus cells with fungal elicitors and analyzed the transcriptome by RNA sequencing. We selected TFs that were significantly increased upon fungal elicitor treatment, and then examined their gene expression by qPCR, and also performed protein network analysis. The results indicated that one unigene (c133966_g3_i1) belonging to the MYB family might play as a master regulator involved in pinosylvin stilbene synthesis under fungal elicitor treatment. In the future, the transgenic pine plants by overexpression of this gene are expected to increase pinosylvin stilbene synthesis which may confer resistance to pine wilt disease.

Appendix

Fig. 12

Annotated unigene ratio (percent) using searched programs of KEGG, NT, EggNOG, Pfam, NR, UniProt, and GO using BLASTN of NCBI BLAST and BLASTX of DIAMOND software with an E-value default cutoff of 1.0E − 5

Fig. 13

Data transformation and normalization with edgeR R library and differentially expressed gene analysis results. a Boxplot of expression difference among six samples. Figures were presented as raw data (Count), Log 2 (counts per Million reads + 1), and TMM Normalization values. A visual representation of a box plot with percentiles (median, 25%, 75%, maximum, and minimum values). b density plot of expression difference among six samples. Raw data (Count), Log 2 (counts per million reads + 1), TMM Normalization values. c Number of significant contigs by fold change and p value. d Log2-fold change and p value obtained from the comparison between two groups plotted as a volcano plot. e A smear plot to show an expression difference in cells between the control and fungal elicitor-treated groups according to the overall average expression level (Smear). (x-axis: average log CPM, y-axis: log2 fold change)

Fig. 14

The PCR efficiency for each primer pair was calculated as E = 10^{(−1/slope)−1}. The logarithm of each concentration of RNA in the dilution series (x-axis) was plotted against the mean of the Ct value for the used concentration (y-axis). E qPCR efficiency, R² linear correlation coefficient

Table 2 Overall mapping ratio for each sample in transcriptome sequencing data of P. strobus cells treated with or without fungal elicitor

Table 3 List of primers for qRT-PCR analysis of P. strobus cells