Comparative Transcriptome Analyses Reveal Different Regulatory Mechanisms in Ecological Adaptation between Chrysanthemum vestitum and Chrysanthemum mongolicum

: Chrysanthemum mongolicum and Chrysanthemum vestitum belong to the Asteraceae family, which occupy a wider range of ecological niches and represent extensive biotic resistance and abiotic tolerance. However, the genetic information about these two species is poorly known, which restricts their utilization. Here, the leaf transcriptomes of the two Chrysanthemum species were investigated and compared. In total, 21,856 differentially expressed genes (DEGs) were identiﬁed between C. vestitum and C. mongolicum , of which 12,076 DEGs were up-regulated and 9780 were down regulated in C. vestitum compared to C. mongolicum . Functional enrichment analysis revealed that stress resistance categories had signiﬁcant proportions. The up-regulated DEGs related to “ABC trans-porters”, “Flavonoid biosynthesis” and “Monoterpenoid biosynthesis” were signiﬁcantly enriched in C. vestitum compared with C. mongolicum . While the DEGs involved in “Biosynthesis of unsaturated fatty acids”, “Proteasome”, “Phenylpropanoid biosynthesis”, “Oxidative phosphorylation”, “Plant-pathogen interaction”, “Starch and sucrose metabolism”, “Glutathione metabolism”, as well as “MAPK signaling pathway” were mostly up-regulated in C. mongolicum compared with C. vestitum , suggesting their important roles in C. mongolicum . These results might explain the differences in morphology and provide potential molecular mechanisms for the ecological adaptation of the two Chrysanthemum species in extreme environments. Together, the results of this study provide a genetic resource that may greatly beneﬁt the genetic improvement of cultivated chrysanthemums and will be helpful for plant conservation and sustainable utilization in the future.


Introduction
Chrysanthemum L. belongs to the Asteraceae family, which comprises about 37 species. More than 20 wild species of Chrysanthemum are native to China, distributed from the north to the south and occurring at elevations ranging from plain to high-alpine, even more than 3600 m, employing a broader diversity of morphological and physiological adaptations to environmental conditions [1,2], and harboring valuable gene resources. Chrysanthemums are one of the most famous floral plants in China and comprise the most abundant cultivars worldwide. However, in the long-term cultivation process, the ornamental characteristics, adaptability, and stress resistance of chrysanthemum gradually decline, which limits its garden application and commercial production. Therefore, the improvement of chrysanthemums with excellent environmental tolerance has always been a hot issue in modern chrysanthemum breeding [3]. Chrysanthemums have been suggested as a complex hybrid that originated from some wild relatives [4,5]. Thus, wild Chrysanthemum species are the primary gene pool of chrysanthemums. The genetic Chrysanthemum species are the primary gene pool of chrysanthemums. The genetic improvement of cultivated chrysanthemums will most benefit from the valuable genetic resources of wild species. However, to date, only a small number of wild Chrysanthemum species genomes have been explored [6][7][8].
Chrysanthemum mongolicum and Chrysanthemum vestitum are two Chrysanthemum wild species that are distributed in arid regions, and have notable tolerance to abiotic stress and barren soil [9,10]. In addition, C. vestitum is a wild hexaploid, which has been supported as wild progenitors that contributed to cultivated chrysanthemums [4,5]. However, the two species morphological traits are clearly distinguishable (Figures 1A, G). Previous studies on these two plants mainly focused on phylogenetic relationships [2,5,10] or on geographical distribution analysis [9,11]. Little is known about the morphophysiological and molecular biology of the overall plant or the plant s resistance characteristics against stress.
. The transcriptome is an efficient method for the specific characterization and quantification of different types of RNA molecules in a tissue of a given organism using highthroughput technologies [12]. It can generate abundant transcripts associated with trait evolution and diversification [13,14], which enable us to understand gene regulation The transcriptome is an efficient method for the specific characterization and quantification of different types of RNA molecules in a tissue of a given organism using highthroughput technologies [12]. It can generate abundant transcripts associated with trait evolution and diversification [13,14], which enable us to understand gene regulation mechanisms and networks and discover interesting genes. In addition, comparative transcriptome studies between closely related species are also important to discover genes related to particular biological functions and provide genus-specific genomic resources. The transcriptome has been extensively used in plants, especially for some non-model species and some species with large and complex genomes, greatly prompting the discovery of novel genes and understanding the complex regulation networks in higher plants [15][16][17].
Here, we characterized the C. mongolicum and C. vestitum transcriptomes using the Illumina paired-end sequencing platform and performed a comparative transcriptome analysis between these two species. Our aims are to find the genetic difference and detect the mechanisms contributing to the strategies they adapt to the environmental stress between C. mongolicum and C. vestitum. The transcriptomic data yielded in the present study provides new information about the adaptation of the two Chrysanthemum species in extreme environments and will be useful for further studies on the transcriptomics and metabolomics of Chrysanthemum and Asteraceae in general. It will also increase the genetic resources available for chrysanthemum breeding or evolutionary analysis.

Plant Materials
C. mongolicum (M) and C. vestitum (V) were collected from Baotou in Inner Mongolia, China, and Funiu Mountain in Henan province, China, respectively. The plants were planted in arid soil filled with stones and bricks in the nurse garden of Northeastern University in Shenyang, China. The voucher specimens were stored in the herbariums of Northeastern University and Northwest A&F University (WUK, Yangling, Shaanxi, China). Fresh leaves of two species were sampled from the transplants being cultivated for one year, immediately frozen in liquid nitrogen, and stored at −80 • C until further use. Three biological replications were created for each species.

Morphological Analysis
The mature leaves of C. mongolicum and C. vestitum were collected and fixed in FAA (formalin-acetic-ethanol-water = 10:5:50:35) for scanning electron microscopy (SEM) analysis. Firstly, the samples were dissected and dehydrated in an ethanol and isoamyl acetate series. Then, they were treated with critical-point drying in CO 2 and sputter-coated with gold according to the description of Usma [18]. The micrographs were taken using a Hitachi S-3500 scanning electron microscope (Hitachi, Tokyo, Japan). The photographs of mature leaves were taken using a Nikon D7100 digital camera.

RNA Extraction, Library Construction, and Illumina Sequencing
Total RNA was extracted using the Omega Plant RNA Kit (Omega, Norcross, GA, USA) according to the manufacturer's instructions and treated with RNase-free DNase I (Omega) to remove the genomic DNA. A NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and an Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) were used to check the quality of RNA. The mRNA purification, fragmentation, cDNA synthesis, adapter addition, fragment size selection, and PCR amplification were conducted according to the method described by Zhang et al. [7]. The cDNA libraries were sequenced on the Illumina sequencing platform (Illumina HiSeq™ 2500) by Novogene Co., Ltd. (Tianjin, China). The raw data were deposited in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the accession number "Bioproject accession: BioProject ID: PRJNA843897 and PRJNA987857".

Sequence Assemble and DEG Analysis
High-quality clean reads were generated after removing the reads having adaptors, containing an 'N' rate greater than 10%, and low-quality reads containing greater than 50% of low-quality (Q-value ≤ 20) bases using the Filterfq program (BGI, Shenzhen, China). The genome of Chrysanthemum nankingense [6] was used as a reference. The clean data were mapped to the C. nankingense genome using Hisat2 software [19]. The differentially expressed genes (DEGs) between the two species were screened using the DESeq R package (version 1.10.1) [20] with a threshold of |log2 (fold change)| ≥ 1 and a false discovery rate (FDR) < 0.05. The accumulation patterns of the DEGs were performed by the heatmap package of the R language. GO and KEGG pathway enrichment analyses of the DEGs were performed with the GOseq R package [21] and KOBAS software [22], respectively.

Validation of DEG Expression
Twelve DEGs were randomly selected to validate the accuracy of RNA-seq data by qRT-PCR. Total RNA was extracted using the E.Z.N.A. Plant RNA Kit (Omega Bio-Tek, Norcross, GA, USA) and digested using RNase-free DNase I (Omega Bio-Tek, Norcross, GA, USA) to remove genomic DNA contamination. Single-strand cDNA was synthesized using a Revert Aid First Strand cDNA Synthesis kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) following the manufacturer's instructions. Specific primers designed by Primer Express 3.0 are listed in Table S1. qRT-PCR was performed in triplicate, containing a template of 1: 5 dilution of the cDNA, 1.0 µL of each primer to a final concentration of 400 nM, 10 µL of SYBR Green Master Mix (Thermo Fisher Scientific, USA), and PCR-grade ddH 2 O to a final volume of 20 µL. The qPCR conditions were carried out as follows: denaturation at 55 • C for 2 min, 94 • C for 5 min, and 40 cycles consisting of 94 • C for 15 s and 60 • C for 1 min on an ABI QuantStudio 5 real-time PCR system thermocycler. Actin was used as an internal control for normalizing the relative expression levels of selected genes according to previously published studies [23,24]. The relative expression levels of DEGs were determined using the 2 −∆∆Ct method [25].

SEM Analysis of Leaf Surface Micro-Morphology of C. mongolicum and C. vestitum
SEM was performed to examine the leaf epidermis micromorphology differences between C. mongolicum and C. vestitum ( Figure 1). The leaf surfaces of the two species showed a distinct difference in appearance. The thick trichome layer is a significant feature of C. vestitum leaves. Both the abaxial and adaxial surfaces of the leaves were covered with densely and thickly appressed pubescence in C. vestitum ( Figure 1A,D,E), in contrast, sparsely short trichomes were present on the leaf surface of C. mongolicum. The stomatal apparatuses and cell shapes were different between the two species. Anomocytic stomata were present in C. vestitum, while C. mongolicum showed actinocytic stomata. The stomata in C. vestitum were larger than those in C. mongolicum, while they were slightly sunken in C. mongolicum ( Figure 1B,C,F,I,L). Leaf trichomes act as a physical barrier, which plays an important role in defending against biotic attack, reducing leaf transpiration, and resisting the harmful effects of UV-B radiation [26][27][28][29]. The stress-tolerant species always present dense trichomes on their leaves. In addition, wax has been suggested to increase the drought tolerance of plants [30,31]. The marked cuticular wax was distinct on the surface of the leaves of C. mongolicum. Thus, these leaf features provide C. vestitum and C. mongolicum with specific eco-physiological advantages to cope with stress. The difference in morphology between C. vestitum and C. mongolicum also revealed that they have evolved distinct adaptation strategies to the conditions of their habitats.

High-Throughput Transcriptome Sequencing and Read Assembly
The leaf transcriptomes of two Chrysanthemum species were characterized using the RNA sequencing method. Approximately 44.43 Gb of Illumina RNA-seq data were obtained from the leaves of C. mongolicum and C. vestitum after the terminal adaptor sequences, duplicated and ambiguous sequences, and low-quality reads were removed. In total, 23.07 Gb and 24.36 Gb of clean read data in the transcriptome were assemblies for C. mongolicum and C. vestitum, respectively. The Q20 scores were greater than 97.7%, and the GC content ranged from 42.69% to 43.1% for each sample (Table 1). A total of 36,333 transcripts were obtained after clean reads of the two species were assembled to the C. nankingense genome sequence, of which 21,129 transcripts were shared between C. vestitum and C. mongolicum, and 8602 and 6602 were specially expressed in C. vestitum and C. mongolicum, respectively ( Figure S1).

Differential Gene Expression Analysis
By filtering with FDR < 0.05 and fold change ≥ 1, 21,856 DEGs in C. vestitum had changed significantly compared with C. mongolicum. Hierarchical clustering analysis (HCA) showed a distinct difference in gene expression patterns between C. vestitum and C. mongolicum ( Figure 2A). There were 12,076 DEGs that were up-regulated in C. vestitum and 9780 that were down-regulated ( Figure 2B). total, 23.07 Gb and 24.36 Gb of clean read data in the transcriptome were assemblies for C. mongolicum and C. vestitum, respectively. The Q20 scores were greater than 97.7%, and the GC content ranged from 42.69% to 43.1% for each sample (Table 1).
A total of 36,333 transcripts were obtained after clean reads of the two species were assembled to the C. nankingense genome sequence, of which 21,129 transcripts were shared between C. vestitum and C. mongolicum, and 8602 and 6602 were specially expressed in C. vestitum and C. mongolicum, respectively ( Figure S1).

Differential Gene Expression Analysis
By filtering with FDR < 0.05 and fold change ≥ 1, 21,856 DEGs in C. vestitum had changed significantly compared with C. mongolicum. Hierarchical clustering analysis (HCA) showed a distinct difference in gene expression patterns between C. vestitum and C. mongolicum (Figure 2A). There were 12,076 DEGs that were up-regulated in C. vestitum and 9780 that were down-regulated ( Figure 2B). .

Enrichment Analysis of DEG
Based on the Nr annotation, the DEGs were classified into 1201 GO categories and mainly assigned to molecular functions and biological processes. In the biological process category, the largest identified term was "response to stress" (GO: 0006950), which contains 216 genes, followed by the aminoglycan catabolic process, amino sugar metabolic

Enrichment Analysis of DEG
Based on the Nr annotation, the DEGs were classified into 1201 GO categories and mainly assigned to molecular functions and biological processes. In the biological process category, the largest identified term was "response to stress" (GO: 0006950), which contains 216 genes, followed by the aminoglycan catabolic process, amino sugar metabolic process, chitin metabolic process, amino sugar catabolic process, and so on. ADP binding, oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular oxygen to two molecules of water, hydrolase activity, hydrolyzing O-glycosyl compounds (GO: 0004553), and hydrolase activity, acting on glycosyl bonds (GO: 0016798), were the most prominent terms in molecular function. Only "nucleolus" was enrichment in the cellular component category ( Figure 3A).
Horticulturae 2023, 9, x FOR PEER REVIEW 6 of 14 process, chitin metabolic process, amino sugar catabolic process, and so on. ADP binding, oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular oxygen to two molecules of water, hydrolase activity, hydrolyzing O-glycosyl compounds (GO: 0004553), and hydrolase activity, acting on glycosyl bonds (GO: 0016798), were the most prominent terms in molecular function. Only "nucleolus" was enrichment in the cellular component category ( Figure 3A). A total of 137 KEGG pathways were enriched among the DEGs between the two species. Among the top 20 enrichment pathways, "Starch and sucrose metabolism" (114), "ABC transporter (111), "Phenylpropanoid biosynthesis (104), "Plant-pathogen interaction"(102), " Proteasome" (98), "Oxidative phosphorylation" (92), "Flavonoid biosynthesis," "Monoterpenoid biosynthesis," "typtophan metabolism," "Zeatin biosynthesis," "Glutathione metabolism," "Cutin suberine and wax biosynthesis," and "Fatty acid metabolism" were most enriched, which were also well-known pathways related to abiotic stress in plants ( Figure 3B). Of which, the up-regulated DEGs related to "ABC transporters", "Flavonoid biosynthesis" and "Monoterpenoid biosynthesis" were significantly enriched in C. vestitum compared to C. mongolicum. While the down-regulated DEGs were significantly enriched in "Biosynthesis of unsaturated fatty acids," "Proteasome," "Phenylpropanoid biosynthesis," "Plant-pathogen interaction, "Oxidative phosphorylation," "MAPK signaling pathway," as well as "Glutathione metabolism" (Figure S2), which indicate the important roles they played in C. mongolicum. These results indicated that different metabolites were synthesized in various Chrysanthemum wild species, which contribute to their morphology and adaptation strategies to their habitats.

Gene Validation and Expression Analysis
Twelve DEGs were randomly chosen for validation of the reliability of transcriptome data using qRT-PCR with gene-specific primers (Table S1). The results of qRT-PCR were significantly consistent with the FPKM value of the transcriptome data, which demonstrated the reliability of the RNA-seq data (Figure 4). A total of 137 KEGG pathways were enriched among the DEGs between the two species. Among the top 20 enrichment pathways, "Starch and sucrose metabolism" (114), "ABC transporter" (111), "Phenylpropanoid biosynthesis" (104), "Plant-pathogen interaction" (102), "Proteasome" (98), "Oxidative phosphorylation" (92), "Flavonoid biosynthesis", "Monoterpenoid biosynthesis", "typtophan metabolism", "Zeatin biosynthesis", "Glutathione metabolism", "Cutin suberine and wax biosynthesis", and "Fatty acid metabolism" were most enriched, which were also well-known pathways related to abiotic stress in plants ( Figure 3B). Of which, the up-regulated DEGs related to "ABC transporters", "Flavonoid biosynthesis" and "Monoterpenoid biosynthesis" were significantly enriched in C. vestitum compared to C. mongolicum. While the down-regulated DEGs were significantly enriched in "Biosynthesis of unsaturated fatty acids", "Proteasome", "Phenylpropanoid biosynthesis", "Plant-pathogen interaction, "Oxidative phosphorylation", "MAPK signaling pathway", as well as "Glutathione metabolism" (Figure S2), which indicate the important roles they played in C. mongolicum. These results indicated that different metabolites were synthesized in various Chrysanthemum wild species, which contribute to their morphology and adaptation strategies to their habitats.

Gene Validation and Expression Analysis
Twelve DEGs were randomly chosen for validation of the reliability of transcriptome data using qRT-PCR with gene-specific primers (Table S1). The results of qRT-PCR were significantly consistent with the FPKM value of the transcriptome data, which demonstrated the reliability of the RNA-seq data (Figure 4).

DEGs Involved in Flavonoid Biosynthesis
Flavonoids are the major secondary metabolites that contribute to plant development, color, and biotic and abiotic stress resistance in plants [31][32][33][34][35]. In this study, a total of 25 DEGs encoding 8 enzymes were assigned to the flavonoid biosynthetic pathway based on a KEGG pathway enrichment analysis. Compared with C. mongolicum, 18 DEGs, including the genes encoding chalcone synthase (CHS), chalcone-flavonone isomerase (CHI), flavonoid 3 -monooxygenase (CYP75B1), and flavone synthase (CYP93B2), were significantly upregulated in C. vestitum. while naringenin 3-dioxygenase (F3H), dihydroflavonol 4-reductase (DFR), and flavonol synthase (FLS) were significantly upregulated in C. mongolicum ( Figure 5, Table S2), which is consistent with the result of Duan et al. [36] that flavonoids were highly accumulated in the leaves of C. mongolicum by metabolism analysis. Flavonoid biosynthesis-associated genes have shown a differential expression pattern in other Chrysanthemum species in previous studies [37]. These results suggested the diversity of flavonoids in Chrysanthemum and that flavonoids might play an important role in the Chrysanthemum species ability to adapt to different environmental conditions.

DEGs Involved in Flavonoid Biosynthesis
Flavonoids are the major secondary metabolites that contribute to plant development, color, and biotic and abiotic stress resistance in plants [31][32][33][34][35]. In this study, a total of 25 DEGs encoding 8 enzymes were assigned to the flavonoid biosynthetic pathway based on a KEGG pathway enrichment analysis. Compared with C. mongolicum, 18 DEGs, including the genes encoding chalcone synthase (CHS), chalcone-flavonone isomerase (CHI), flavonoid 3'-monooxygenase (CYP75B1), and flavone synthase (CYP93B2), were significantly upregulated in C. vestitum. while naringenin 3-dioxygenase (F3H), dihydroflavonol 4-reductase (DFR), and flavonol synthase (FLS) were significantly upregulated in C. mongolicum ( Figure 5, Table S2), which is consistent with the result of Duan et al. [36] that flavonoids were highly accumulated in the leaves of C. mongolicum by metabolism analysis. Flavonoid biosynthesis-associated genes have shown a differential expression pattern in other Chrysanthemum species in previous studies [37]. These results suggested the diversity of flavonoids in Chrysanthemum and that flavonoids might play an important role in the Chrysanthemum species' ability to adapt to different environmental conditions. culturae 2023, 9, x FOR PEER REVIEW 8 of Figure 5. Expression patterns of genes involved in the flavonoid biosynthesis pathway. The expr sion levels of DEGs in the heatmap obtained with Log2FC. Red represents up−regulated genes, a blue represents down−regulated genes. One the left was C. vestitum and on the right was C. mong icum in heat map. KEGG Pathway was derived from https://www.kegg.jp/kegg/kegg1.html ( cessed on 25 April 2023).

Candidate genes involved in terpenoid biosynthesis
Terpenoids are the most structurally and quantitatively abundant compounds pr duced by plants [38], which are the major compounds of aromatics and play an importa role in plant defense against environmental stresses [39]. In the present study, the u regulated DEGs related to "Monoterpenoid biosynthesis" were significantly enriched C. vestitum compared to C. mongolicum. All of the DEGs assigned to encoding saluta dine reductase (SalR), neomenthol dehydrogenase (NDS), and 10-hydroxygeraniol deh drogenase (10HGO), which were key genes related to monoterpene biosynthetic, we upregulated in C. vestitum (Table S3). Trichomes have been reported to be able to synth size and accumulate terpene compounds [40][41][42]. The highly expressed genes related monoterpene biosynthesis were consistent with the densely covered trichomes in t leaves of C. vestitum, suggesting large amounts of bioactive monoterpenes in its leav The DEGs related to terpenoid biosynthesis have been reported to be higher expressed C. indicum var. aromaticum than in C. nankingense and C. indicum [37]. These results ind cated that terpenoid content varies in Chrysanthemum species and that C. mongolicu might be a potential resource for the aromatic improvement of chrysanthemums.

Candidate Genes Involved in Terpenoid Biosynthesis
Terpenoids are the most structurally and quantitatively abundant compounds produced by plants [38], which are the major compounds of aromatics and play an important role in plant defense against environmental stresses [39]. In the present study, the upregulated DEGs related to "Monoterpenoid biosynthesis" were significantly enriched in C. vestitum compared to C. mongolicum. All of the DEGs assigned to encoding salutaridine reductase (SalR), neomenthol dehydrogenase (NDS), and 10-hydroxygeraniol dehydrogenase (10HGO), which were key genes related to monoterpene biosynthetic, were upregulated in C. vestitum (Table S3). Trichomes have been reported to be able to synthesize and accumulate terpene compounds [40][41][42]. The highly expressed genes related to monoterpene biosynthesis were consistent with the densely covered trichomes in the leaves of C. vestitum, suggesting large amounts of bioactive monoterpenes in its leaves. The DEGs related to terpenoid biosynthesis have been reported to be higher expressed in C. indicum var. aromaticum than in C. nankingense and C. indicum [37]. These results indicated that terpenoid content varies in Chrysanthemum species and that C. mongolicum might be a potential resource for the aromatic improvement of chrysanthemums.

Candidate Genes Involved in Starch and Sucrose Metabolism Pathway
Starch and sucrose metabolism provide fuel for plant growth and development [43,44]. The accumulation of soluble sugars is also an important stress signal for plants [45,46]. Here, the starch and sucrose metabolism pathways were significantly enriched, and 114 DEGs involved in these pathways were detected (Table S4), of which 61 DEGs, including the genes encoding granule-bound starch synthase I (WAXY), starch synthase (glgA), glucose-1phosphate adenylyltransferase (glgC), Amylose (AMY), amylomaltase (malQ), and sucrose synthase (SUS), were upregulated in C. vestitum, while 58 were up-regulated in C. mongolicum, including the genes encoding trehalose-6-phosphate synthase (TPS), Sucrose Phosphatase (SPP), Trehalase (TREH), etc. The results indicated that the two species contain abundant but different soluble sugars that help them adapt to the harsh environment.

Biosynthesis of Unsaturated Fatty Acids
Unsaturated fatty acids (UFAs) as membrane ingredients and carbon and energy resources in triacylglycerols (TAGs) are indispensable compounds for higher organisms [47,48]. In addition, C18 UFAs function as antioxidants and barrier constituents and play crucial roles in plant defense against both abiotic and biotic stresses [49,50]. Fatty acid desaturase-2 (FAD2) is a key gene converting oleic acid (18:1) to linoleic acid (18:2) [51,52]. In our study, the biosynthesis of the unsaturated fatty acids pathway was significantly enriched, and 38 DEGs were detected, of which 32 DEGs encoding fatty acid desaturase-related genes were upregulated in C. mongolicum (Table S5), suggesting unsaturated fatty acids may be responsible for the resistance to abiotic stress in this species.

Candidate Genes Involved in ABC Transporters
The ATP-binding cassette (ABC) transporters, which served as exporters and importers, are involved in regulating plant development, nutrient uptake, osmotic homeostasis maintenance, pathogen resistance, metal tolerance, etc. [53][54][55]. The ABC transporter pathway was significantly enriched, and 111 DEGs were identified in our study, of which 70 DEGs were up-regulated in C. vestitum, including PDR5, ABCA, ABCB, ABCC, and ABCG (Table S6). The ABCG proteins have been suggested as candidates for heavy metal transporters, e.g., cadmium, lead, etc. [56,57], suggesting C. vestitum might be a potential plant for restoration of heavy metal-contamination.

DEGs Involved in Plant-Pathogen Interaction
A multi-layered innate immune system has been constructed in higher-land plants to reject or attenuate pathogen infection. PAMP-triggered immunity (PTI) and effectortriggered immunity (ETI) are two major innate immune systems in response to pathogen stress [58]. They can induce transcriptional reprogramming of different sets of defenserelated genes to confer resistance. Calcium signaling has been reported to respond to various biotic and abiotic stresses in plants [59]. In the present study, 102 DEGs involved in the plant-pathogen interaction pathway were identified (Table S7). Thirty-nine DEGs, including PTI6, PTI1, and WRKY, showed higher expression in C. vestitum than in C. mongolicum, while 63 DEGs, including 22 genes encoding a calcium-binding protein (CML), 2 genes encoding cyclic nucleotide gated channel (CNGC), 18 DEGs encoding calciumdependent protein kinase (CPK), two RPS2, and SGT1, were significantly upregulated in C. mongolicum. These results indicated that the PTI reaction and calcium channel played important roles in the plant-pathogen interactions of C. vestitum and C. mongolicum, respectively.

DEGs Involved in Phenylpropanoid Biosynthesis
Phenylpropane compounds are widely distributed in plants. They play an important role in plant development, flower pigments, and helping plants defend against various environmental stresses [34,60]. In the present study, 104 DEGs involved in the Phenylpropanoid biosynthesis pathway were detected ( Figure 6; Table S8) between C. vestitum and C. mongolicum. Most of the DEGs (64) in the phenylpropanoid biosynthesis pathway were up-regulated in C. mongolicum, including the genes encoding cinnamyl-alcohol dehydrogenase (CAD), caffeoyl-CoA O-methyltransferase (CCOMT), Peroxidase (PER), shikimate O-hydroxycinnamoyl-transferase (HST), coniferyl-aldehyde dehydrogenase (REF), coniferyl-alcohol glucosyltransferase (UGT72E), ferulate-5-hydroxylase (CYP84A), caffeoyl shikimateesterase (CSE), and cinnamyl-alcohol dehydrogenase that related to ligin biosynthesis (Table S8). Lignin has important roles in plant resistance against various biotic and abiotic stresses [61]. This result indicated that high lignification occurred in the leaves of C. mongolicum to adapt to the adverse environmental conditions, which agrees with the appearance that the leaves at the branch bottom were withered and shrunk in this species. lturae 2023, 9, x FOR PEER REVIEW 10 shikimateesterase (CSE), and cinnamyl-alcohol dehydrogenase that related to ligin synthesis (Table S8). Lignin has important roles in plant resistance against various b and abiotic stresses [61]. This result indicated that high lignification occurred in the le of C. mongolicum to adapt to the adverse environmental conditions, which agrees with appearance that the leaves at the branch bottom were withered and shrunk in this spe

Transcription Factors Analysis of C. vestitum and C. mongolicum
Transcription factors (TFs) play a crucial role in regulating plant development environmental stress responses [62]. Numerous stress-related TFs have been identifie various plants. In this study, 802 differentially expressed TFs belonging to 78 TF fam were identified. AP2/ERF is the most predominant category with 63 members identi followed by MYB (57), bHLH (49), NAC (47), C2H2 (43), and WRKY (35) (Table  AP2/ERF, MYB, NAC, bHLH, and WRKY are notable TFs involved in the regulatio responses to abiotic stress in plants [63][64][65][66][67]. The special function of these TFs in these species need to be further investigated with a functional genomics approach.

Conclusions
In the present study, the leaf transcriptomes of the two wild Chrysanthemum sp were compared, and a large number of genomic resources for chrysanthemums wer tained. A distinctly differentiated gene expression pattern was found between C. vest and C. mongolicum, and the differences and complexities of secondary metabolic bio

Transcription Factors Analysis of C. vestitum and C. mongolicum
Transcription factors (TFs) play a crucial role in regulating plant development and environmental stress responses [62]. Numerous stress-related TFs have been identified in various plants. In this study, 802 differentially expressed TFs belonging to 78 TF families were identified. AP2/ERF is the most predominant category with 63 members identified, followed by MYB (57), bHLH (49), NAC (47), C2H2 (43), and WRKY (35) (Table S9). AP2/ERF, MYB, NAC, bHLH, and WRKY are notable TFs involved in the regulation of responses to abiotic stress in plants [63][64][65][66][67]. The special function of these TFs in these two species need to be further investigated with a functional genomics approach.

Conclusions
In the present study, the leaf transcriptomes of the two wild Chrysanthemum species were compared, and a large number of genomic resources for chrysanthemums were obtained. A distinctly differentiated gene expression pattern was found between C. vestitum and C. mongolicum, and the differences and complexities of secondary metabolic biosynthesis were discovered in the two Chrysanthemum species. The DEGs associated with "ABC transporters", "Flavonoid biosynthesis" and "Monoterpenoid biosynthesis" might play a crucial role in adaption to environmental stresses for C. vestitum, while the genes related to "Biosynthesis of unsaturated fatty acids", "Proteasome", "Phenylpropanoid biosynthesis", "Oxidative phosphorylation", "Plant-pathogen interaction", "Starch and sucrose metabolism", "Glutathione metabolism", as well as "MAPK signaling pathway" might be responsible for adaption of the adverse condition in C. mongolicum. Thus, C. mongolicum and C. vestitum have different stress resistance and positive adaptation mechanisms. The large EST sequences generated in this study help us better understand the genetic characteristics and resistance mechanisms of wild chrysanthemum plants. However, identifying conserved and divergent genes involved in ecological adaptation will be helpful to understand the evolutionary processes underlying the adaptation of C. vestitum and C. mongolicum. Moreover, with the increasing demand for chrysanthemums, there is an urgent need to develop attractive new varieties with extensive resistance to different environmental stresses.
The key candidate genes involved in various pathways that respond to stress found here need to be further explored by transgenic or gene-editing techniques, which will accelerate the genetic improvement of cultivated chrysanthemums.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/horticulturae9080868/s1. Table S1: Primers of differential expressed genes used for qRT-PCR; Table S2: Genes involved in Flavonoid biosynthesis; Table S3: Genes involved in terpenoid biosynthesis; Table S4: DEGs related to starch and sucrose metabolism pathway; Table S5: Biosynthesis of unsaturated fatty acids; Table S6: Genes involved in ABC transporters; Table S7: DEGs involved in plant-pathogen interaction; Table S8: DEGs involved in phenylpropanoid biosynthesis; Table S9: DEGs encoding transcription factors analysis between C. vestitum and C. mongolicum; Figure S1: Venn diagram of special and common transcripts between C. vestitum and C. mongolicum. CV: C. vestitum; CM: C. mongolicum; Figure S2: KEGG enrichment analysis of top 20 differentially expressed genes between C. vestitum and C. mongolicum. A represents the enriched pathways of up-regulated genes and B is the enriched pathways of down-regulated genes. The dot size means the number of DEGs in a related pathway. The rich factor represents the ratio of the DEG numbers assigned to a certain pathway to all gene numbers annotated in this pathway.

Data Availability Statement:
The data that support the findings of this study are openly available in the GenBank of NCBI at https://www.ncbi.nlm.nih.gov/. The accession numbers can be found in the article.

Conflicts of Interest:
There are no conflicts of interest exist in the submission of this manuscript.