The Trichoderma longibrachiatum MD33 was isolated from stem segments of a single wild D. nobile (Sarsaiya et al., 2019a,2020b) and stored at the laboratory of Nanjing Tech University’s Institute of Bioresources Engineering. After activating the stored fungus for 48 h, it was transferred to PD medium and cultivated at 28°C for 36 h. Thereafter 0.1180 g MeJA (95%) was completely dissolved in pure ethanol, transferred to a 10 ml volumetric flask and dilute with deionized water to 50 mmol/L MeJA. Methyl Jasmonate (MeJA) was added to PD medium at 0, 5, 20, and 35 μmol/L (250 ml Erlenmeyer flask, 100 ml PD medium), eight bottles per group. The pH value at the beginning was neutral. The block at the colony’s edge where the new hypha forms was removed and put into each flask using a 6 mm punch. The culture conditions used with the strain were neutral pH, 28°C, 120 rpm rotation speed, and 15 days (Sarsaiya et al., 2020a). The growth of MD33 were shown in Figure 1. There were three biological replicates for each group.

Five bottles of Trichoderma longibrachiatum MD33 were used for dendrobine extraction and detection. The following steps were used: The fungi cultured in shake flasks were filtered, dried to constant weight at 45°C, and ground to a powder. Thereafter 50 ml chloroform was added to the powder, transferred to a separation funnel for static stratification after 30 min of ultrasonic extraction, distilled and the lower solvent layer was dried under reduced pressure (rotary evaporator SHZ-III model). The dried powder was dissolved in 5 ml methanol centrifuged at 1,000 rpm for 15 min. Thereafter 1 ml solution was filtered via a 0.22 μm filter membrane, and a sample was injected for detection.

The LC-MS technique was used to characterize non-volatile and thermally fragile compounds. The detection of dendrobine was performed using a UHPLC system (Thermo Fisher Scientific Dionex Ultimate 3000, Golden Valley, Minnesota, United States) with a column (150 × 2.1 mm,1.9 μm) and a mobile phase consisting of 0.1% formic acid: acetonitrile at 95:5 (v/v) with a flow rate of 0.3 ml/min, column temperature of 40°C, feed volume 2 μl, a sheath gas flow rate of 35 arbitrary units, an auxiliary gas flow rate of 15 arbitrary units, a spray voltage of 3.5 kV, a capillary temperature of 350°C, an aux gas heater temperature of 300°C, the detection range was 100–1,500 m/Z, the scanning mode was full scan/positive ion scan, and the resolution was MS full scan 70,000 full width at half maxima (FWHM).

To provide a scientific rational for MeJA-induced dendrobine accumulation in MD33, a standard curve was established and the contents of dendrobine were evaluated under control and MeJA conditions. For the standard curve, 1.00 mg of dendrobine standard was accurately weighed, dissolved in methanol to prepare 1 g/L mother liquor, and then diluted into standard solutions at concentrations of 1, 10, 100, 500, 1,000, and 1,500 ng/L. The linear regression equation and standard curve were constructed using the LC-MS detection results. The regression curve equation of dendrobine is y = 0.6726x + 23.205, R² = 0.9948 (Figure 2A). Each sample detected for three times.

The cultivated Trichoderma longibrachiatum MD33 was treated according to the LC/MS findings (Step 2.1 and Step 2.2). Three flasks from each group were randomly selected and filtered. Prior to RNA extraction, the hypha samples were snap-frozen in liquid nitrogen and stored at –80°C.

Total RNA was extracted from hypha using the RNeasy mini kit (QIAGEN, Cat. No. 74104, Germany) and genomic DNA was removed using DNase. Then RNA quality was determined by Bioanalyzer 2100 and the RNA 6000 Nano LabChip kit (Agilent, CA, United States). Only high-quality RNA sample (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5) was used A SureSelect strand-specific RNA library preparation kit from Agilent Technologies was used to generate an RNA-seq library according to the manufacturer’s procedure using 2 μg of total RNA. The Agilent 2100 Bioanalyzer and the ABI Step One Plus Real-Time PCR system was used to assess the quality of the RNA-seq library. Thereafter, sequencing was performed on the Illumina HiSeq 4000 platform. The transcriptome library was prepared by KEGENE Company (Shandong, China). There were three biological replicates for each group.

The raw FASTQ data was analyzed for quality using FastQC and filtered to remove reads with sequencing adaptors, unknown nucleotides (Nradio > 5%), and poor quality (quality scores 30). Due to the unavailability of a reference genomic sequence, de novo assembly was used to create the transcripts. The redundant sequences were eliminated using the CORSET software suite, and the longest unigenes were further spliced. All assembled unigenes were aligned against the non-redundant (Nr), nucleotide (Nt), protein family (Pfam), Gene Ontol-ogy (GO), SwissProt, and KOG/COG, Kyoto Encyclopedia of Genes and Genomes (KEGG) databases with a threshold of E-value < 10^–10.

In this study, the expression level of all transcripts was standardized using RSEM and Bowtie2 with default parameters. Ballgown was used for Differentially Expressed Gene (DEG) investigation. An FDR value of (0.05) and a log2FC value of 1 were used as the criteria for DEG screening. The DEGs were then analyzed for GO and KEGG enrichment.

qRT-PCR was performed to validate the accuracy of RNA-seq data. GAPDH gene was used as an internal reference. The foldchange in gene expression was calculated using the comparative Ct method (1 + E^–△△Ct) (Pfaffl, 2001). Total RNA samples were extracted from the hyphae using RNeasy mini kit (QIAGEN, Cat: 74104, Germany). The cDNA synthesis kit (TSINGKE, TSK302M) was used to synthesize cDNA from 1 μg of RNA. For gene expression analysis, the Applied Biosystems StepOnePlus™ Real-Time PCR system was used. qRT-PCR amplification was performed in 20 μl reactions containing 10 μl 2 × TSINKE^® Master qPCR Mix (SYBR GREEN), 20 ng cDNA, and 0.4 μM of each primer (Supplementary Table 5). The PCR reaction conditions were: 95°C for 15 min, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s. After final annealing (72°C, 5 min) and redenaturation (95°C, 30 s), a melt curve analysis was done by increasing from 60 to 95°C at 0.5°C intervals.

The LC-MS results revealed that standard chemical reference dendrobine was found at a retention time of 3.05 and a molecular weight of 264.195. From the MD33 dendrobine the dendrobine peak was recorded at 3.07 with a molecular weight of 264.195 (Figures 2B,C). LC-MS detection results showed that, the content of dendrobine in MeJA treatment groups were higher than control group. When the concentration of MeJA was 20 μmol/L, the content of dendrobine was the highest (0.68 μg/L), which increased 44.7% compared with the control group (0.47 μg/L) (Figure 3A). indicating that MeJA plays an important role in dendrobine accumulation in MD33.

Three distinct RNA-seq libraries were generated utilizing RNA from the control and MeJA treated groups. A total of 274,332,098 pure reads were isolated from 288,958,158 raw reads in this study, equal to 41.16 Gb of clean bases (Supplementary Table 1). The Q20/Q30 fraction, the N proportion, and the GC percentages were 98.73/95.64, 0.00, and 58.05%, respectively. The clean sequences were then de novo assembled and clustered into 21,331 (median length = 3,850 bp) and 4,477 (median length = 3,326 bp) genes, respectively, using Trinity (Supplementary Tables 2, 3). We analyzed the size distributions of transcripts and genes, respectively. For transcripts, the majority (76.22%) were > 2,000 bp in length, 23.00% were 500–2,000 bp in length, and 0.77% were 500 bp in length. For genes, 42.93% were > 2,000 bp in length, 54.75% were 500–2,000 bp on length, and 2.32% < 500 bp in length (Supplementary Table 2). All identified genes were annotated using several databases. The NR database contained 6,037 genes, the KO database contained 1,765 genes, the KEGG database contained 1,761 genes, the KO database contained 2,920 genes, the SwissProt database contained 4,885 genes, the Pfam database contained 5,666 genes, the GO database contained 4,677 genes, and the KOG database contained 5,777 genes (Supplementary Table 4).

To ascertain the likely function of these DEGs in Trichoderma longibrachiatum MD33, we performed a Gene Ontol-ogy (GO) enrichment analysis, which revealed that these DEGs were classified into three primary categories: biological process (BP), cellular component (CC), and molecular function (MF). The terms “oxidation-reduction process,” “protein transport,” and “cell cycle” were the most frequently used in the biological process category; “nucleus,” “membrane,” and “cytoplasm” were frequently used in the cellular component category; and a large percentage of genes were associated with “metal ion binding,” “nucleotide binding,” and “hydrolase activity” in the molecular function category (Figure 4A). Moreover 1,761 genes were classified into 18 Kyoto Encyclopedia of Genes and Genome (KEGG) secondary pathways. The majority of pathways were classified as “Metabolism” (866 genes) or “Genetic Information Processing” (695 genes). Moreover, 190 genes belonged to the “amino acid metabolism” pathway, 139 genes belonged to the “carbohydrate metabolism” pathway, 126 genes belonged to the “metabolism of cofactors and vitamins” pathway, 106 genes belonged to the “lipid metabolism” pathway, and 28 genes belonged to the “metabolism of terpenoids and polyketides” pathway in metabolism, 302 genes belonged to the “translation” route in “Genetic Information Processing,” whereas 205 genes belonged to the “folding, sorting, and degradation” pathway (Figure 4B).

The DEGs that increased considerably after MeJA treatment were regarded as MeJA-regulated genes (Supplementary Table 6). We performed a significance analysis on the DEGs and displayed the results using a volcano diagram. Notably, 444 DEGs were identified between the control and MeJA treatment groups. Among these significant DEGs, 226 (50.9%) were up-regulated, while the remaining 218 (49.1%) were down-regulated (Figure 3B). The examination of GO enrichment revealed that the most enriched GO keywords were “oxidoreductase activity,” “ion transport,” “multi-organism,” “interspecies interaction between organisms,” “pathogenesis,” “anion transport,” and “ligase activity” (Figure 4C). Additionally, genes were allocated to several KEGG metabolic pathways (Supplementary Table 4). Enrichment analysis showed that the significantly enriched KEGG pathways were “Translation,” “Amino acid metabolism,” “Folding,” “sorting and degradation,” “Carbohydrate metabolism,” “Transport and catabolism,” “Lipid metabolism,” “Metabolism of cofactors and vitamins” (Figure 4D).

The DEGs analysis of secondary metabolites revealed that, 10 clusters were up-regulated compares to the control group. These clusters are involved in secondary metabolite production and degradation, caffeine metabolism, glucosinolate biosynthesis, and novobiocin biosynthesis. There was a significant difference in Phenylpropanoid production (p < 0.05), and the function of the divergent gene may be beta-glucosidase. The unchanged clusters were involved in the synthesis of tropane, piperidine, and pyridine alkaloid (Supplementary Table 5).

The upstream biosynthesis routes for sesquiterpene intermediate products have been extensively researched, and are preserved in plants (Li et al., 2017). Dendrobium alkaloids are mainly sesquiterpene alkaloids, which are derived from MVA, and the MEP route would act as a source of isoprene units in the production of dendrobine (Chen et al., 2019). Three DEGs, isopentenyl-diphosphate Delta isomerase (IDI), dishomevalonate decarboxylase (MVD), and farnesyl diphosphate synthase (FDPS), were mapped to the MVA pathway in this study (Figure 5).

By examining the DEGs, we identified 51 potential P450 unigenes, nine of which were novel. MD33 with and without MeJA treatment is shown in Table 1. Cluster-4787.0, Cluster-514.0, Cluster-2347, Cluster-5538.0, and Cluster-5724.0 were found to be up-regulated, while Cluster-1305.0, Cluster-3418.0, Cluster-1102.0, and Cluster-3538.0 were found to be down-regulated (Figure 6 and Table 1).

Aminotransferase and methyltransferase were investigated as two additional key enzymes involved in post-modification. Clustering analysis indicated that only one aminotransferase (Cluster-6002.0) was down-regulated (p < 0.05), but Cluster-4588.0 (log2 fold change = 3.16-fold), Cluster-5638.0 (log2 fold change = 1.85-fold), and Cluster-5199.0 (log2 fold change = 3.02-fold) were up-regulated after MD33 MeJA treatment (Table 1). According to the reference genome (T. longibrachiatum ATCC 18648) (Shi et al., 2020), Cluster-6002.0 encodes a glutamate-1-semialdehyde aminotransferase. The metabolic pathways implicated may include heme (k01845) production and porphyrin and chlorophyll metabolism (ko00860). Heme and porphyrin are pyrrole chemicals with a similar structure to dendrobine. Hence this gene may be implicated in the transaminase function during dendrobine production (Table 1).

By analyzing the DEGs, we obtained 140 putative methyltransferase unigenes, of which eight unigenes (Cluster-1958.0, Cluster-2098.0, and Cluster-4527.0 increased expression, Cluster-3608.0, Cluster-3923.0, Cluster-4296.0, Cluster-2665.0, and Cluster-3774.0 decreased expression) expressed differently (p < 0.05) between MeJA treatment MD33 and control MD33 (Figure 6 and Supplementary Table 1). The primary molecular functions were hydrolase and methyltransferase activities.

Various Transcription Factors (TFs) were reported to participate in the biosynthesis of secondary metabolites in plants. In this study, 484 putative TF genes belonging to 31 major TF families were identified, and 17 genes belonging to 11 major TF families were analyzed (Figure 6). The Zn-clus (six genes) and the zinc finger protein contained the maximum members. Among these TF genes, three Zn-Clus genes, one SET gene, one C3H gene, and 1 GNAT gene were found to be up-regulated in the MeJA treatment MD33 group (Table 1).

A qRT-PCR assay with independent samples from the control and MeJA treatment groups was used to verify the expression changes of several key genes involved in the dendrobine biosynthesis pathway. In total 10 genes, including one MVA pathway genes, two P450 family genes, two methyltransferase genes, three aminotransferase genes, and two TF genes, were selected to confirm the RNA-seq data. The expression levels of these selected genes were consistent with the RNA-seq results (Figure 7 and Supplementary Table 5).