Comparative Transcriptomic Analysis of Genes in the 20-Hydroxyecdysone Biosynthesis in the Fern Microsorum scolopendria towards Challenges with Foliar Application of Chitosan

Microsorum scolopendria is an important medicinal plant that belongs to the Polypodiaceae family. In this study, we analyzed the effects of foliar spraying of chitosan on growth promotion and 20-hydroxyecdysone (20E) production in M. scolopendria. Treatment with chitosan at a concentration of 50 mg/L in both young and mature sterile fronds induced the highest increase in the amount of accumulated 20E. Using RNA sequencing, we identified 3552 differentially expressed genes (DEGs) in response to chitosan treatment. The identified DEGs were associated with 236 metabolic pathways. We identified several DEGs involved in the terpenoid and steroid biosynthetic pathways that might be associated with secondary metabolite 20E biosynthesis. Eight upregulated genes involved in cholesterol and phytosterol biosynthetic pathway, five upregulated genes related to the methylerythritol 4-phosphate (MEP) and mevalonate (MVA) pathways, and several DEGs that are members of cytochrome P450s and ABC transporters were identified. Quantitative real-time RT-PCR confirmed the results of RNA-sequencing. Taken together, we showed that chitosan treatment increased plant dry weight and 20E accumulation in M. scolopendria. RNA-sequencing and DEG analyses revealed key enzymes that might be related to the production of the secondary metabolite 20E in M. scolopendria.


Introduction
Microsorum scolopendria commonly known as "wart fern" is a valuable air-purifying ornamental plant and an important medicinal plant belonging to the Polypodiaceae family [1]. In traditional Polynesian medicine, M. scolopendria fronds are used as antioxidants and as medication for numerous physical and mental disorders [2,3]. The adaptogenic and anabolic effects of M. scolopendria demonstrated by modern pharmacological research are attributed to 20-hydroxyecdysone (20E), a non-androgenic steroid used as an ingredient in dietary supplements [4,5] to increase energy and muscle mass [6]. M. scolopendria fronds and 20E exert a neuroprotective effect in aging rats [2,6], improve cognitive impairment [7], protect against brain injury, and inhibit the production of reactive oxygen species In the present study, we used chitosan as an elicitor to reveal the molecular effects of chitosan in stimulating the synthesis of secondary metabolites in M. scolopendria fronds [25]. The fronds of M. scolopendria were sprayed with chitosan at three different concentrations, to test the effect of chitosan on plant growth promotion. After chitosan treatment, foliar application of chitosan positively affected the growth of M. scolopendria by increasing the dried mass of the fronds. However, the fern fronds treated with chitosan had slightly lower fresh weight than that of the control, whereas the chitosan-treated plants showed increased dry weight 7-14 days after treatment (DAT) compared with that of the control (Figure 1). The plants sprayed with chitosan at 50 mg/mL exhibited the highest dry weight in both young (1.93 times) and mature (1.21 times) sterile fronds after 7 DAT compared with that of the control (Figure 1). Specifically, fern plants sprayed with chitosan at 50 and 100 mg/mL significantly enhanced the growth rate of young sterile fronds at 7 DAT (by 12.3% and 11.4%, respectively) when compared with that at 0 DAT, whereas the growth rate of mature sterile fronds slightly increased by 5.6% and 4.8%, respectively. However, the effect of 200 mg/mL chitosan application did not significantly increase the relative growth rate (RGR), which was only slightly enhanced by 0.84 times in the mature sterile fronds when compared with that of the control (Figure 1). These data indicated that plants treated with chitosan at concentrations of 50 and 100 mg/mL showed increased growth compared with that in the control based on dried mass accumulation. Of the three different chitosan concentrations, chitosan treatment at 50 mg/mL showed the highest value for dried mass. We used 50 mg/mL as the optimal concentration for further experiments.

20E Content Increased in Sterile Fronds of Chitosan-treated M. scolopendria
Regarding plant growth and its attributes, chitosan elicits the activation of enzymes in the metabolic pathway linked with numerous secondary metabolites. Our previous study showed 20E accumulation in different frond tissues of the fern Microsorum; however, sterile fronds showed the highest 20E accumulation [26]. Thus, we focused on the application of chitosan to elucidate 20E and the associated mechanisms affecting plant growth and 20E biosynthesis. Our results showed that the plants treated with 50 mg/mL chitosan had significantly higher 20E production in young and mature sterile fronds compared with that in pre-treated fronds. As shown in Figure 2, the highest 20E production (59.14 mg/g) was observed in the mature sterile fronds at 14 DAT. However, the 20E content in the young and mature sterile fronds increased 5.79 and 3.60 times at 7 DAT and 8.11 and 3.87 times at 14 DAT, respectively, compared with that on 0 DAT ( Figure 2). Therefore, foliar spraying of 50 mg/mL chitosan enhanced the 20E content in both young and mature sterile fronds.

20E Content Increased in Sterile Fronds of Chitosan-Treated M. scolopendria
Regarding plant growth and its attributes, chitosan elicits the activation of enzymes in the metabolic pathway linked with numerous secondary metabolites. Our previous study showed 20E accumulation in different frond tissues of the fern Microsorum; however, sterile fronds showed the highest 20E accumulation [26]. Thus, we focused on the application of chitosan to elucidate 20E and the associated mechanisms affecting plant growth and 20E biosynthesis. Our results showed that the plants treated with 50 mg/mL chitosan had significantly higher 20E production in young and mature sterile fronds compared with that in pre-treated fronds. As shown in Figure 2, the highest 20E production (59.14 mg/g) was observed in the mature sterile fronds at 14 DAT. However, the 20E content in the young and mature sterile fronds increased 5.79 and 3.60 times at 7 DAT and 8.11 and 3.87 times at 14 DAT, respectively, compared with that on 0 DAT ( Figure 2). Therefore, foliar spraying of 50 mg/mL chitosan enhanced the 20E content in both young and mature sterile fronds.

Identification of Differential Expressed Genes (DEGs) of M. scolopendria in Response to Chitosan Treatment
To narrow down significant DEGs, we applied fold change > 2 and adjusted p-value < 0.01 as cutoffs resulting in 121 DEGs that were further divided into 61 upregulated and 60 downregulated genes (Table S1 and Figure 3A). For instance, representative highly upregulated genes were those encoding DMR6-LIKE OXYGENASE 2-like protein, MADSbox protein, DMR6-LIKE OXYGENASE 2-like protein, L-type lectin-domain containing receptor kinase IX.1-like protein, and terpene synthase 4 (Table S1).
The highly downregulated genes were those encoding lectin 2a, low molecular mass early light-inducible protein HV90, COBRA-like protein, MFT-like protein, reticulon-like protein B10, C2 domain-containing protein, and PAR1-like protein (Table S1). Hierarchical clustering using the 121 DEGs was conducted in six different samples and shows two distinct groups of samples and two groups of genes ( Figure 3B).

Identification of Differential Expressed Genes (DEGs) of M. scolopendria in Response to Chitosan Treatment
To narrow down significant DEGs, we applied fold change > 2 and adjusted p-value < 0.01 as cutoffs resulting in 121 DEGs that were further divided into 61 upregulated and 60 downregulated genes (Table S1 and Figure 3A). For instance, representative highly upregulated genes were those encoding DMR6-LIKE OXYGENASE 2-like protein, MADS-box protein, DMR6-LIKE OXYGENASE 2-like protein, L-type lectin-domain containing receptor kinase IX.1-like protein, and terpene synthase 4 (Table S1).
The highly downregulated genes were those encoding lectin 2a, low molecular mass early light-inducible protein HV90, COBRA-like protein, MFT-like protein, reticulon-like protein B10, C2 domain-containing protein, and PAR1-like protein (Table S1). Hierarchical clustering using the 121 DEGs was conducted in six different samples and shows two distinct groups of samples and two groups of genes ( Figure 3B).

Functional Annotation and Classification of Assembled Transcripts
To further determine the enriched functions of the identified 3552 DEGs, we performed a gene ontology (GO) enrichment analysis using Fisher's exact test with an FDRadjusted p-value < 0.01 as the cutoff. All DEGs were functionally categorized based on GO terms according to three functional categories: biological process, cellular component, and molecular function. Significantly enriched GO terms were visualized (Table S2 and Figure  4A). We identified 18 enriched GO terms according to cellular components. Many DEGs were associated with a cellular anatomical entity (59 DEGs), intracellular (27 DEGs), membrane (13 DEGs), an integral component of membrane (9 DEGs), an intrinsic component of membrane (9 DEGs), which were found to be specific to secondary metabolite production. According to molecular function, 29 GO terms were enriched. Of these, GO terms associated with catalytic activity (29 DEGs), binding (21 DEGs), transporter activity (7 DEGs), and transferase activity (7 DEGs) were significantly enriched (Table S3 and Figure  4A). According to the biological processes, 139 GO terms were enriched and associated with secondary metabolites and plant defense mechanisms. For instance, the significantly enriched GO terms were associated with metabolic processes (32 DEGs), followed by cellular processes (27 DEGs), triterpenoid biosynthetic processes (eight DEGs), organic substance metabolic processes (six DEGs), cellular metabolic processes (five DEGs), primary metabolic processes (five DEGs), and oxidation-reduction processes (four DEGs). In addition, GO terms associated with responses to lipid and carbohydrate metabolic processes were significantly enriched ( Figure 4A).
To obtain a more detailed insight into biological function, GO enrichment analysis was performed on a dataset of significant 121 DEGs. Using a threshold of false discovery rate (FDR)<0.05, we found that many DEGs were targeted to three cellular components:

Functional Annotation and Classification of Assembled Transcripts
To further determine the enriched functions of the identified 3552 DEGs, we performed a gene ontology (GO) enrichment analysis using Fisher's exact test with an FDR-adjusted p-value < 0.01 as the cutoff. All DEGs were functionally categorized based on GO terms according to three functional categories: biological process, cellular component, and molecular function. Significantly enriched GO terms were visualized (Table S2 and Figure 4A). We identified 18 enriched GO terms according to cellular components. Many DEGs were associated with a cellular anatomical entity (59 DEGs), intracellular (27 DEGs), membrane (13 DEGs), an integral component of membrane (9 DEGs), an intrinsic component of membrane (9 DEGs), which were found to be specific to secondary metabolite production. According to molecular function, 29 GO terms were enriched. Of these, GO terms associated with catalytic activity (29 DEGs), binding (21 DEGs), transporter activity (7 DEGs), and transferase activity (7 DEGs) were significantly enriched (Table S3 and Figure 4A). According to the biological processes, 139 GO terms were enriched and associated with secondary metabolites and plant defense mechanisms. For instance, the significantly enriched GO terms were associated with metabolic processes (32 DEGs), followed by cellular processes (27 DEGs), triterpenoid biosynthetic processes (eight DEGs), organic substance metabolic processes (six DEGs), cellular metabolic processes (five DEGs), primary metabolic processes (five DEGs), and oxidation-reduction processes (four DEGs). In addition, GO terms associated with responses to lipid and carbohydrate metabolic processes were significantly enriched ( Figure 4A). ribosomes, membranes, and integral components of the membrane. According to the biological process, the top five enriched GO terms were associated with transmembrane transport, translation, response to bacteria, biosynthetic processes, and protein phosphorylation. The enriched GO terms according to molecular function were ATP binding, oxidoreductase activity, iron ion binding, heme binding, and monooxygenase activity.  To obtain a more detailed insight into biological function, GO enrichment analysis was performed on a dataset of significant 121 DEGs. Using a threshold of false discovery rate (FDR) < 0.05, we found that many DEGs were targeted to three cellular components: ribosomes, membranes, and integral components of the membrane. According to the biological process, the top five enriched GO terms were associated with transmembrane transport, translation, response to bacteria, biosynthetic processes, and protein phosphorylation. The enriched GO terms according to molecular function were ATP binding, oxidoreductase activity, iron ion binding, heme binding, and monooxygenase activity.
Next, we performed an enrichment analysis against the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway to identify the key metabolic pathways enriched in the DEGs. All identified DEGs were assigned to at least one KEGG pathway and categorized into 236 metabolic pathways (Table S3). Of these, 38 pathways such as ABC transport, starch and sucrose metabolism, metabolic pathway, amino sugar and nucleotide sugar metabolism, and MAPK signaling pathway were significantly overrepresented under chitosan treatment (Table S3).
We also found that several DEG-encoding enzymes are involved in the steroid and phytosterol biosynthetic pathways for 20E synthesis ( Figure 6). For example, we identified several DEGs encoding eight enzymes in the steroid biosynthesis pathway, mainly enriched in the cholesterol and phytosterol biosynthetic pathways. Moreover, we identified several DEGs encoding sterol-4alpha-carboxylate 3-dehydrogenase (E1. 1    The functional domains of the proteins encoded by the upregulated genes were classified according to the InterPro database. As a result, a total of 41 protein domains were significantly enriched. The most enriched functional categories were transporter, signal transduction, and cytochrome P450 domains ( Figure 7A). In contrast, only 17 protein domains were identified in downregulated genes. They were reticulon, stay-green protein, and plant lipid transfer protein domains ( Figure 7B). Additionally, several transcription factors belonging to the AP2/ERF, MYB, and WRKY families were identified, which are involved in plant development and lipid metabolism (Table S3 and Figure 7). Many DEGs were also associated with transmembrane transport and ATP binding ( Figure S1). Notably, numerous genes (51 DEGs) encoding cytochrome P450s were identified (Table S3). Based on functional annotation, various cytochrome P450s (cytochrome P450 71, 72, 83, and 709) related to specialized triterpene metabolism were identified ( Table 2 and Table S3). All the DEGs encoding cytochrome P450s were up-regulated ( Figure S1). The expression levels of CYP72A15 (2.164), CYP709B2 (2.227), CYP83B1 (4.452), and CYP71A1 (5.679) increased after chitosan treatment (Table 1). The seven DEGs including cytochrome P450 genes (CYP72A15 and CYP709B2), triterpenoid biosynthetic genes (ADH and MADS-box), ABC transporter-type G gene, and biosynthetic process-related genes (reticulon and LTP) were selected for the validation of RNA-seq results using quantitative real-time RT-PCR. We found that the expression tendency of all examined genes was consistent between the RNA-seq and qRT-PCR results. For example, genes encoding ABC-type G, ADH, MADS-box transcription factor, CYP72A15, and CYP709B2 were upregulated, whereas two genes encoding reticulon and LTP were downregulated ( Figure 8).   biosynthetic process-related genes (reticulon and LTP) were selected for the validation RNA-seq results using quantitative real-time RT-PCR. We found that the expression te dency of all examined genes was consistent between the RNA-seq and qRT-PCR resul For example, genes encoding ABC-type G, ADH, MADS-box transcription facto CYP72A15, and CYP709B2 were upregulated, whereas two genes encoding reticulon an LTP were downregulated (Figure 8).

Discussion
The detectable amounts of 20E in plants are very low, accounting for 1-2% of the d weight in the few identified high-accumulators, whereas most 20E-accumulating speci contain only 0.01-0.1% [27]. The 20E content must be considered not only when selectin a high-producing genetic cultivar, but also when considering treating the plants with a propriate elicitors. In this study, we examined the effects of foliar spraying of chitosan o the physio-biochemical and molecular responses of plants. Our results showed that a three concentrations (50, 100, and 200 mg/mL) of chitosan significantly (p<0.05) affecte the growth of M. scolopendria. Of the three different concentrations of chitosan, 50 mg/m showed the highest increase in RGR in young and mature sterile fronds at 7 DAT. Mor over, chitosan treatment provided the highest growth-stimulating effect in young steri fronds compared with that of the control. Natural polysaccharides, such as chitosan an its derivatives, improved the growth attributes and efficiency of uptake of individual el ments, which in turn is associated with the intensification of plant growth [28,29]. Sever

Discussion
The detectable amounts of 20E in plants are very low, accounting for 1-2% of the dry weight in the few identified high-accumulators, whereas most 20E-accumulating species contain only 0.01-0.1% [27]. The 20E content must be considered not only when selecting a high-producing genetic cultivar, but also when considering treating the plants with appropriate elicitors. In this study, we examined the effects of foliar spraying of chitosan on the physio-biochemical and molecular responses of plants. Our results showed that all three concentrations (50, 100, and 200 mg/mL) of chitosan significantly (p < 0.05) affected the growth of M. scolopendria. Of the three different concentrations of chitosan, 50 mg/mL showed the highest increase in RGR in young and mature sterile fronds at 7 DAT. Moreover, chitosan treatment provided the highest growth-stimulating effect in young sterile fronds compared with that of the control. Natural polysaccharides, such as chitosan and its derivatives, improved the growth attributes and efficiency of uptake of individual elements, which in turn is associated with the intensification of plant growth [28,29]. Several previous studies have shown that chitosan treatment at lower concentrations results in increased vegetative growth in marigold, grape, okra, freesia, basil, and strawberry [30][31][32][33][34][35]. Moreover, the foliar application of chitosan inhibited the transpiration and gaseous exchange and increased the absorption of water and principal nutrients, thus providing water to biomass and yield. For example, the foliar application of chitosan solution to Dracocephalum kotschyi resulted in a significant increase in biomass [18,36]. Our study confirmed previous findings and validated that the application of chitosan promotes plant growth. Furthermore, our results demonstrated that elicitor concentration and timing of application might be very important in determining the optimal effects of chitosan.
Interestingly, foliar application of chitosan improved 20E accumulation in the sterile fronds of M. scolopendria compared with that in the untreated control. Chitosan foliar spray in common rue plants led to the accumulation of various types of secondary metabolites such as coumarins, acridone, quinolone alkaloids, and flavonoids [37]. Chitosan application is an effective elicitor for improving rosmarinic acid, quercetin, and apigenin levels up to 16-fold higher than that of the control [18]. Many investigators have reported the properties of chitosan in plant growth induction and secondary metabolite elicitation. However, only a few studies provide limited data on the effects of spray application of chitosan in the fern M. scolopendria, which remarkably increased the active 20E content in sterile fronds. This study shows that low doses of chitosan (50 mg/mL) sprayed on the fronds of M. scolopendria increased growth with a concurrent significant enhancement of 20E content in sterile fronds.
Transcriptome analysis of mature sterile fronds of M. scolopendria has been carried out as the fronds are the primary site of synthesis and accumulation of bioactive 20E. Chitosan is the main signal of secondary metabolite production across the plant kingdom, from angiosperms to tracheophytes. Chitosan treatment triggers the biosynthesis of the majority of secondary metabolites (i.e., terpenoids and phenylpropanoids) through extensive transcriptional reprogramming. The elicitation of secondary metabolite 20E enriched in M. scolopendria might be due to the role of chitosan in gene regulation and signaling, which leads to the activation of enzymes in metabolic pathways linked with 20E biosynthesis. In this study, using RNA-seq, we investigated transcriptional changes in M. scolopendria in response to 50 mg/mL chitosan. As a result, we narrowed the number of DEGs from 3552 to 121 ones. Similarly, several significant DEGs have been identified in response to chitosan in potato leaves, suggesting that the activation of a set of conserved genes positively influences plant growth, development, and metabolism [38]. Chitosan treatment induced the expression of 83 DEGs, most of which are involved in biological processes that promote photosynthesis and respiration. When we compared the RNA-seq data of the chitosan-treated ferns and potatoes, transcripts encoding mitochondrial proteins were shown to distinguish between these two treated plants, such as the genes encoding cytochrome c oxidase and ATP synthase subunits, which were upregulated in potatoes [38], whereas cytochrome c oxidase was not influenced by chitosan, and ATP synthase subunits were downregulated in M. scolopendria.
The identification of candidate genes and key enzymes is crucial for understanding the biosynthetic pathways of active 20E in M. scolopendria. 20E is a terpenoid synthesized from the terpenoid backbone biosynthesis pathway, which includes both MEP and MVA pathways [39]. The precursor of 20E biosynthesis operates more readily and intracellularly by following the plastidial MEP pathway and cytosolic MVA pathway, respectively [10]. Based on the enriched GO terms according to cellular components, both the plasticidal MEP and cytosolic MVA pathways were functional for the biosynthesis of 20E in Microsorum. Moreover, analyses of genes against the KEGG pathway revealed that many of the identified DEGs were involved in 236 metabolic pathways. They encode various enzymes and proteins involved in the biosynthesis of different isoprenoids such as monoterpenes, diterpenes, triterpenes, sesquiterpenoids, and steroids.
A previous report also predicted that the major role of the terpenoid backbone biosynthesis pathway is specialized secondary metabolite biosynthesis [40]. The MVA pathway generally supplies the precursors for the production of sesquiterpenes, triterpenes, and brassinosteroids. The MEP pathway generally supplies the precursors for the biosynthesis of diterpenoids [41]. Identification of genes encoding 3-hydroxy-3-methylglutaryl coenzyme A reductase (HMGR) in this study correlated with an earlier study demonstrating that this enzyme is the rate-limiting enzyme of the MVA and cholesterol pathways. Overall, four HMGR genes were identified in the mature sterile fronds of M. scolopendria after foliar chitosan treatment. We noticed the involvement of HMGR in MVA. Moreover, cholesterol induced by chitosan treatment was required for the upstream processes of 20E biosynthesis. Additionally, we identified transferases, such as 2-C-methyl-Derythritol-4-phosphate cytidylyltransferase (ispD), sterol 24-C-methyltransferase (ERG6), 24-methylenesterol C-methyltransferase, and reductases, such as 7-dehydrocholesterol reductase (DHCR) and hydroxymethylglutaryl-CoA reductase (HMGR). Genes associated with the triterpenoid biosynthesis pathway have frequently been identified in several medicinal plants [42]. Moreover, genes encoding key enzymes involved upstream of the 20E biosynthetic pathway, such as steroid 17α-monooxygenase, farnesyl diphosphate synthase, farnesol kinase, prenylcysteine α-carboxymethylesterase, diphosphomevalonate decarboxylase, 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase, and 2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase were identified in the present study.
We identified 64 genes encoding ABC transporters. Genes in the ABC transporter family, based on the hydrolysis of ATP, are known to play a major role in the transport of secondary metabolites [43,44]. Since many plant secondary metabolites are used medicinally, ABC transporters are also associated with drug resistance [45]. For example, the indirect role of ABC transporters in regulating secondary metabolism has been studied by overexpressing enzymes involved in the primary metabolism [46]. Meanwhile, an inspection of the leaf epidermis enriched transcription database identified an ABC transporter that mediates the transport of anticancer drug components in the Catharanthus [45]. Therefore, ABC transporters may be associated with the transport of special compounds, including alkaloids, terpenoids, and polyphenols [47], and the secondary metabolite 20E.
We identified several transcription factors involved in plant development and lipid metabolism including AP2/ERF, MYB, and WRKY. The ethylene responsive factor (AP2/ERF) superfamily and R2R3-MYB family belong to one of the largest diverse families of transcription factors and are involved in controlling plant lipid metabolism processes, biosynthesis of primary and secondary metabolites, and adaptation to stress [48]. Previous studies have shown that WRKY and ERF transcription factors regulate sterol synthesis in plants. WRKY1 belongs to the WRKY TFs and participates in regulating the levels of secondary metabolites such as sterols and glycoalkaloids in plants. For example, WsWRKY1 silencing resulted in the hindrance of plant growth and decreased plant sterol content, whereas the overexpression of WsWRKY1 promoted the synthesis of plant sterols in tobacco and tomato [49]. Furthermore, with the overexpression of PqWRKY1 in Arabidopsis, the expression of plant sterol synthesis-related genes (such as HMGR, FPS2, SQS1, and SQE2) in transgenic lines increased one-to five-fold compared with that in the control [50].
Several upregulated genes encode amino acid transporters, which are essential for many substrates in plants and are involved in shuttling many secondary metabolite intermediates between subcellular compartments. The cytochrome P450 superfamily mediates the oxidation of various substrates involved in both primary and secondary metabolism, which enhances 20E accumulation in chitosan-treated M. scolopendria [51,52]. Recently, 52 genes of the cytochrome P450s family have been identified, which play crucial roles in terpenoid biosynthesis in various medicinal plants. The majority of these belong to the CYP71 family [53]. The most downregulated genes encoded members of the reticulon-like protein family and glutamate dehydrogenase (Table S4). In plants, reticulons are integral ER membrane proteins that are often associated with lipid synthesis, which might have a direct negative influence on 20E biosynthesis in M. scolopendria [54].
Based on gene annotation analysis of all DEGs, low concentrations (50 mg/mL) of chitosan application act as plant signal molecules to improve plant growth by inducing secondary metabolite 20E accumulation in the medicinal fern M. scolopendria. This process might occur via several mechanisms, such as (1) triggering various important transportermediated membrane transport, such as ATP-binding cassette (ABC) transporters and major facilitator superfamily (MFS) transporters, and repression of lipid transfer proteins that play a role in lipid transfer across the cytoplasm in plants (we suppose that these transporters play important roles linked to chitosan application in 20E biosynthesis by transporting substrate primary metabolites, especially lipid and terpenoid) and (2) stimulating plant steroid biosynthesis via highly expressed cytochrome P450s genes for a role in 20E biosynthesis. A previous study screened the five most highly expressed cytochrome P450 genes for phytoecdysteroid biosynthesis in the hairy roots of Ajuga [55]. The proposed model for the influence of chitosan on 20E elicitation is depicted in Figure 7C.
In conclusion, chitosan treatment increased the dry weight of the plant and 20E accumulation in M. scolopendria. RNA sequencing and DEG analyses revealed key enzymes that might be related to the production of the secondary metabolite 20E in M. scolopendria. We believe that the identified key genes associated with chitosan signaling and 20E production in M. scolopendria could be useful for further genetic analysis of 20E biosynthesis and pathway engineering.

Plant Materials and Foliar Chitosan Application
The fern M. scolopendria (Burm f.) Copel was grown in compost soil without fertilizer and under natural conditions in the experimental plot in a greenhouse at Kasetsart University, KamphaengSaen campus, Nakhon Pathom, Thailand. Each pot was arranged at a density of two pots per square meter. Chitosan from shrimp shells with an average molecular weight of 100 kDa (85% deacetylated chitin) (Sigma-Aldrich, Saint Louis, MO, USA) was dissolved in 5% acetic acid and diluted in distilled water to three different concentrations (50, 100, and 200 mg/mL). The pH of the chitosan solutions was adjusted to 6.5 with 2 M NaOH. Two months after transplantation, the plants were separated into four groups and arranged in a completely randomized block design with three replicates. One untreated group was designated as the control plants, whereas the remaining three groups were sprayed with chitosan solutions of various concentrations using an ultra-low volume sprayer. The spray application of chitosan on fully expanded fronds was performed in the early morning of every alternate day and continued for two weeks. The number of spray shots per plant was maintained at approximately 100 mL per pot. Young and mature sterile fronds were collected on day 0 after treatment and at 7 and 14 DAT. The fronds were rinsed with distilled water and dried at 60 • C until a constant weight was achieved. For transcriptome analyses, mature sterile fronds were sampled, immediately frozen, and stored at -80 • C.

Plant Growth Measurement
Plant growth was measured based on the fresh and dry weights of the fronds. At three different time points (0, 7, and 14 DAT), three different concentrations of chitosan treatment (50, 100, and 200 mg/mL) were applied. Subsequently, young and mature sterile fronds were harvested and weighed before and after drying at 50 • C in a hot-air oven until a constant weight was achieved. Growth rates at consecutive treatment intervals were measured or estimated to determine the RGR. RGR (g g −1 d −1 ) was calculated as (InW 2 − InW 1 )/(t 2 − t 1 ) for each individual M. scolopendria, where W 1 and W 2 are the mean dry mass of individual plants at time t 1 (the initial day of treatment) and time t 2 (the interval between treatments, days), respectively.

Determination of 20-Hydroxyecdysone Content
The M. scolopendria plants were separated into two groups and arranged in a randomized block design with three replicates. One group was the control with no foliar spraying, whereas the other group was sprayed with 50 mg/mL chitosan using the method described earlier. Young and mature sterile fronds were harvested at 0, 7, and 14 d after chitosan treatment and dried in a forced-air oven at 50 • C until a constant weight was achieved. Dry fronds were powdered using a grinder and stored at -20 • C until extraction. The 20-hydroxyecdysone content was determined using a method described previously [1]. Dried tissue powder (200 g) was extracted with 2 L ethanol at 28-30 • C for 48 h. The supernatant was evaporated until the sample was dry, which was re-dissolved in methanol (5 mL). After filtration using a nylon membrane filter (0.45 µm) (Millipore, Darmstadt, Germany), aliquots of the reaction products (20 µL) were injected into a C18 reverse-phase column (4.6 mm × 250 mm) for HPLC analysis (Water Corporation, Miford, MA, USA). The mobile phase was a mixture of acetonitrile:propan-2-ol (5:2, v/v) or trifluoroacetic acid in water (20:80, v/v), at a flow rate of 1 ml min −1 ; UV detection was performed at a wavelength of 245 nm. Data, including the means of three replicates and ± standard error values, are presented as error bars.

Total RNA Extraction, cDNA Library Construction, and RNA Sequencing
RNA-sequencing was conducted to establish the transcriptional profiles of mature sterile fronds after treatment with 50 mg/mL chitosan at 14 DAT. Three replicates from each M. scolopendria frond treatment group were ground into powder in a mortar with liquid nitrogen, and total RNA was extracted from the samples using NucleoSpin ® RNA (MACHEREY-NAGEL, Düren, Germany) according to the manufacturer's protocol. The quality and quantity of RNA were assessed using gel electrophoresis, an ND-1000 Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA), and an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). DNA was removed from the samples using RNase-free DNase (Promega, Madison, WI, USA) according to the manufacturer's instructions [56]. The RNA integrity number (RIN) values of these six samples ranged from eight to nine, which showed good quality for further downstream processes. After that, six cDNA libraries were generated for RNA-seq and paired-end sequenced using the NovaSeq 6000 platform. After removing the adapter sequences, duplicate sequences, ambiguous reads, and low-quality reads, clean reads between 22.18 and 30.74 million paired ends were generated ( Table 2). cDNA synthesis was performed using equal quantities of high-quality RNA from each sample following TruSeq RNA Library Prep Kit v2 (Illumina, San Diego, CA, USA). Sequencing libraries were generated according to the manufacturer's instructions (Illumina, San Diego, CA, USA). For paired-end sequencing, both ends of six cDNA libraries from three replications were sequenced using an Illumina NovaSeq 6000 system (Illumina, San Diego, CA, USA). The quantity and quality of the cDNA library were regulated by Macrogen (Seoul, Republic of Korea).

Assessment of Differential Gene Expression
Raw data obtained from the six samples were deposited in the Sequence Read Archive (SRA) database of NCBI with accession numbers: SRR8480631, SRR8480739, SRR8480744, SRR8481433, SRR8481526, and SRR8481655. High-quality reads were filtered by eliminating low-quality reads (Q < 20) from raw reads using BBDuk (https://jgi.doe.gov/data-andtools/software-tools/bbtools/bb-tools-user-guide/bbduk-guide/ accessed on 1 December 2022). Filtered reads were used for mapping. The transcriptome of M. scolopendria from a previous study was used as a reference transcriptome [26]. Raw data were mapped to the reference transcriptome using the BBMap assembler with default parameters (https://jgi. doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/bbmap-guide/ accessed on 1 December 2022). We obtained fragments per kilobase of transcripts per million mapped reads (FPKM) and mapped the read numbers for each transcript. DEGs were identified by comparing chitosan-treated samples (three samples) to non-treated samples (three samples). The mapped reads were used for DEG analysis using DESeq2 implemented in DEBrowser v1.24.1 [57]. Based on a fold change > two times and adjusted p-value < 0.05, we identified a total of 3552 DEGs in response to chitosan treatment. To identify significant DEGs, we used fold change > two times and adjusted p-value < 0.01 as cutoffs. The expression levels of the individual transcripts were quantified in terms of FPKM values.

Sequencing Annotation and Classification
We used annotated information for M. scolopendria transcriptome from a previous study [26]. For example, we evaluated several databases used previously to annotate M. scolopendria transcripts, including the National Center of Biotechnology Information (NCBI) non-redundant protein (NR, https://www.ncbi.nlm.nih.gov/protein/ accessed on 1 December 2022), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www. genome.jp/kegg/ko.html accessed on 1 December 2022), Gene Ontology (GO, http:// www.geneontology.org/ accessed on 1 December 2022), UniProt (http://www.uniprot.org/ accessed on 1 December 2022), and EggNOG (http://eggnogdb.embl.de/ accessed on 1 December 2022). GO terms were subsequently classified using Blast2GO [58] via OmicsBox with an e-value default cutoff. The enriched GO terms were identified according to three functional categories: biological processes (BP), cellular components (CC), and molecular functions (MF). For KEGG pathway analysis of the DEGs, bidirectional best hit (BBH) was used to search against the KEGG database to obtain the KO (reference pathway) number of KEGG annotations [59].

Validation of the Expression of Selected Genes Using qRT-PCR
We validated the RNA-seq results of seven selected genes involved in phytosterol biosynthesis using quantitative reverse transcription-polymerase chain reaction (qRT-PCR) analysis. Primers were designed using the Primer 3 software (www.embnet.sk/cgi-bin/ primer3_www.cgi accessed on 1 December 2022). The primer names and sequences used for the primer design are listed in Table S5. Total RNA was extracted from each sample using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Total RNA was treated with DNase I (Thermo Scientific, Vilnius, Lithuania), and cDNA synthesis was performed using oligo(dT) primers and the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Dreieich, Germany). qRT-PCR was performed with 1 µL of cDNA in a total volume of 10 µL using iTaq Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) according to the manufacturer's instructions. In detail, 2 µL of cDNA (1:10 diluted) was placed in a 20 µL reaction mixture containing 10 µL of iTaq Universal SYBR Green Supermix (2x) and 0.2 µL of each primer (10 pmol). qRT-PCR was run on a CFX96 qPCR machine (Bio-Rad, Hercules, CA, USA), and the thermal cycling conditions were adjusted as follows: 95 • C for 3 min, followed by 40 cycles of 95 • C for 10 s and 53 • C for 30 s. Actin was used as the reference gene to normalize the gene expression level, and relative gene expression was calculated using the Ct (2 −∆∆Ct ) method [60]. The experiment was performed with three biological replicates. Three technical duplicates were performed for every biological replicate. The outcomes of each biological replicate were calculated using the averages of the three technical replicates.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms24032397/s1. Table S1: The list of differentially expressed genes (DEGs) between chitosan-treated and untreated M. scolopendria at 14 days after treatment (DAT); Table S2: The gene ontology (GO) classification of up-regulated and down-regulated DEGs; Table S3. The Fragments Per Kilobase of transcripts per million mapped reads (FPKM) value and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation of DEGs; Table S4: The list of DEGs encoded for catalytic activity, transmembrane transport, ATP binding and cytochrome P450s; Table S5. Primers for quantitative real-time PCR; Figure S1: Heat map diagram representing transcript expression of DEG function in the most enriched GO categories such as catalytic activity, transmembrane transport, ATP binding, and cytochrome P450s.