Comparative Transcriptome Analysis of Waterlogging-Sensitive and Waterlogging-Tolerant Chrysanthemum morifolium Cultivars under Waterlogging Stress and Reoxygenation Conditions

Waterlogging stress is among the most severe abiotic stressors in the breeding and the production of Chrysanthemum morifolium. However, the mechanism underlying the response to waterlogging and post-waterlogging reoxygenation in C. morifolium remains unknown. In this study, we compared the differences between the transcriptomes of two chrysanthemum cultivars, i.e., the waterlogging-tolerant cultivar “Nannongxuefeng” and the waterlogging-sensitive cultivar “Qinglu”, by performing RNA-seq to elucidate the possible mechanism of waterlogging and reoxygenation in C. morifolium. “Nannongxuefeng” had a higher ethylene production under the waterlogging and reoxygenation conditions. Furthermore, the expression of transcription factors and genes that are involved in the hormone response, N-end rule pathway and ROS signaling significantly differed between the two cultivars. “Nannongxuefeng” and “Qinglu” significantly differed in their response to waterlogging and reoxygenation, providing a deeper understanding of the mechanism underlying the response to waterlogging and guidance for the breeding of C. morifolium.


Introduction
Water is crucial for plant growth and development. However, many plants are often subjected to flooding stress due to excess water remaining in soils. Flooding stress is among the most widespread abiotic stressors during the lifespan of plants [1]. Flooding involves the following two major stressors: submergence, which occurs when the entire plant (root and shoot) is under water, and waterlogging, which occurs when only the root-zone of the plant is flooded [2]. The excess water content in the soil replaces the gas spaces surrounding the roots and shoots, thus rendering the plants hypoxic or anoxic [3]. The low oxygen status of plants caused by flooding results in the suppression of various necessary physiological activities, which jeopardizes plant growth and development [4].
When oxygen is limited, aerobic respiration is suppressed in the roots [4]. Anaerobic respiration becomes the primary pathway supplying energy to the plants. Fermentation pathway enzymes, such as pyruvate decarboxylase (PDC), alcohol dehydrogenase (ADH), and sucrose synthase (SUS), are involved in sugar metabolism and are significantly up-regulated during the waterlogging process [5][6][7][8][9]. The limited Chrysanthemum is a leading ornamental species that is widely popular worldwide [45]. However, most cultivars are susceptible to waterlogging stress and exhibit a wilted, yellow, rotted phenotype under stress conditions. However, a few tolerant cultivars have been identified [46,47]. The morphological, anatomical, physiological, and biochemical responses of chrysanthemum to waterlogging have been described [46,48]. In addition, ethylene production is enhanced by waterlogging in two chrysanthemum cultivars with different tolerance levels, and the tolerant cultivar ("53-4") has a higher ethylene production than the susceptible cultivar ("13-13") [46]. Forty-five molecular markers are significantly associated with the waterlogging tolerance of chrysanthemum, and four favorable parental lines have been identified as potential donors for improving waterlogging tolerance [47]. Nevertheless, the molecular mechanism of the response of chrysanthemum to waterlogging is unclear due to its complicated genetic background and the lack of reference genome information. Here, two chrysanthemum cultivars, i.e., "Qinglu" (waterlogging-sensitive) and "Nannongxuefeng" (waterlogging-tolerant), were examined in order to determine the possible difference in their transcriptome profiles in response to waterlogging stress and post-stress reoxygenation and the signaling pathways and regulatory networks that are related to waterlogging and reoxygenation.

Morphological Changes and Ethylene Measurement in "Nannongxuefeng" and "Qinglu" Following Waterlogging Stress
Here, we compared the morphological changes in two different cultivars, i.e., "Nannongxuefeng" and "Qinglu", following 6 d of waterlogging stress. No morphological changes were observed in "Nannongxuefeng", which displayed green leaves and an upright stem ( Figure  1c,d). However, "Qinglu" displayed a damaged phenotype with a bent stem and wilted leaves (Figure 1a,b). The mature leaves were completely withered and yellow, and only the young leaves remained partially green (Figure 1a,b). Thus, "Nannongxuefeng" was more tolerant than "Qinglu" to waterlogging stress. The ethylene production was measured 0 h, 3 h, 6 h, 12 h, and 24 h after waterlogging and after 2 h of reoxygenation after 12 h of waterlogging of the roots. The following four groups were assessed: control, waterlogging, waterlogging + ACC and waterlogging + AVG. The ethylene production in "Nannongxuefeng" was higher than that in "Qinglu" at each time point in all groups regardless of the waterlogging stage or reoxygenation stage (Table 1). During the waterlogging stage, the waterlogging + ACC group had a higher ethylene level than the other groups, followed by The ethylene production was measured 0 h, 3 h, 6 h, 12 h, and 24 h after waterlogging and after 2 h of reoxygenation after 12 h of waterlogging of the roots. The following four groups were assessed: control, waterlogging, waterlogging + ACC and waterlogging + AVG. The ethylene production in "Nannongxuefeng" was higher than that in "Qinglu" at each time point in all groups regardless of the waterlogging stage or reoxygenation stage (Table 1). During the waterlogging stage, the waterlogging + ACC group had a higher ethylene level than the other groups, followed by the waterlogging group, the control group, and finally, the waterlogging + AVG group at all five time points (Table 1). The same pattern was observed during the reoxygenation stage. Furthermore, the ethylene production peaked at 6 h, declined at 12 h, and dramatically increased at 24 h in both cultivars (Table 1). However, the ethylene production after 12 h + 2 h reoxygenation was higher than that after the 12 h of waterlogging treatment (Table 1). We speculate that the higher ethylene production during waterlogging may be a vital mechanism of the tolerance to waterlogging in "Nannongxuefeng". The increase in the ethylene production after reoxygenation may be due to the ACC accumulation during submergence in the chrysanthemum roots, and the higher oxygen density promoted the transformation of ACC to ethylene [49]. Table 1. Measurement of ethylene production in "Qinglu" and "Nannongxuefeng" at 0 h, 3 h, 6 h, 12 h, 24 h, and 12 h + 2 h recovery. ACC, 1-aminocyclopropane-1-carboxylic; AVG, Aminoethoxyvinyl glycine. Values are presented as the mean ± SE. Letters next to the numbers indicate significant differences between the respective values (p < 0.05).

Gene Annotation and Functional Classification
To compare the differentially expressed genes (DEGs) between the two cultivars in response to waterlogging and reoxygenation, we separated the up-regulated and down-regulated genes in the four groups (QCK 12 h vs. QW (Figure 2A, Tables S1-S4). The number of up-regulated genes in each group was 18,407, 4177, 11,734, and 5128 ( Figure 2A, Tables S1-S4). The number of down-regulated genes in each group was 14,197,53,413,28,745, and 61,689 ( Figure 2A, Tables S1-S4). In these four comparisons, the tolerant cultivar "Nannongxuefeng" had more DEGs than the sensitive cultivar "Qinglu" during both the waterlogging and reoxygenation stages (Figure 2A, Tables S1-S4), demonstrating that "Nannongxuefeng" had a stronger response to waterlogging and reoxygenation than "Qinglu" at the transcription level. The DEGs in the other five comparisons (i.e., Q 0 h vs. X 0 h, QCK 12 h vs. XCK 12 h, QW 12 h vs. XW 12 h, QCK 12 h + 2 h vs. XCK 12 h + 2 h, and QW 12 h + 2 h vs. XW 12 h + 2 h) are also displayed in Figure 2A Venn diagrams of the abovementioned DEGs during the waterlogging and reoxygenation stages were also performed ( Figure 2B,C). In comparison of CK 12 h vs. W 12 h, 6127 up-regulated DEGs were found to be common to both cultivars ( Figure 2B). The sensitive cultivar "Qinglu" showed more number of up-regulated DEGs (12,280 genes) when compared to the tolerant cultivar "Nannongxuefeng" (5607 genes) ( Figure 2B). Among the down-regulated DEGs, 7840 DEGs were found to be common to both the cultivars ( Figure 2B). 6357 DEGs were specifically down-regulated in "Qinglu" and 20,905 were found to be specific to "Nannongxuefeng" ( Figure 2B). In comparison of CK 12 h + 2 h vs. W 12 h + 2 h, 2090 up-regulated DEGs were found to be common to both cultivars ( Figure 2C). 2087 DEGs were specifically up-regulated in "Qinglu" and 3038 were found to be specific to "Nannongxuefeng" ( Figure 2C). Among the down-regulated DEGs, 41,288 DEGs were found to be common to both the cultivars ( Figure 2C), suggesting that down-regulated DEGs overlapped greatly under reoxygenation conditions between two caltivars. However, 12,125 DEGs were specifically down-regulated in "Qinglu", and 20,401 were found to be specific to "Nannongxuefeng" ( Figure 2C). The abovementioned results demonstrated that many DEGs varied between the two cultivars in response to waterlogging and reoxygenation ( Figure 2, Tables S5-S9). Thus, sensitive and tolerant cultivars may have two independent mechanisms for coping with waterlogging stress and short-term recovery. To deeply understand the DEGs between the two cultivars in response to waterlogging stress and reoxygenation, the DEGs were further investigated by performing a GO functional analysis. In the Q 0 h vs. X 0 h comparison, 6921 DEGs were assigned a GO term. Nine GO classes significantly fell into the category 'cellular component', 8 GO classes significantly fell into the category "molecular function", and 8 GO classes significantly fell into the category 'biological process' (Table  S10). In the QCK 12 h vs. XCK 12 h comparison, 16,575 DEGs were assigned a GO term. Twenty-one GO classes significantly fell into the category "cellular component", 46 GO classes significantly fell into the category "molecular function", and 33 GO classes significantly fell into the category "biological process" (Table S11). In the QW 12 h vs. XW 12 h comparison, 7041 DEGs were assigned a GO term. Three GO classes significantly fell into the category "molecular function", and 13 GO classes significantly fell into the category "biological process" (Table S12). In the QCK 12 h + 2 h vs. XCK 12 h + 2 h comparison, 7374 DEGs were assigned a GO term. Fifteen GO classes significantly fell into the category "cellular component", 52 GO classes significantly fell into the category "molecular function", and 34 GO classes significantly fell into the category 'biological process' (Table S13). In the QW 12 h + 2 h vs. XW 12 h + 2 h comparison, 9443 DEGs were assigned a GO term. Seven GO classes significantly fell into the category "cellular component", 30 GO classes significantly fell into the category "molecular function", and 10 GO classes significantly fell into the category "biological process" (Table S14).
A KEGG pathway analysis of the abovementioned DEGs was also performed. In the Q 0 h vs. X 0 h comparison, 7623 DEGs were assigned to 135 pathways in the KEGG database. "DNA replication", "Isoquinoline alkaloid biosynthesis", "Tyrosine metabolism", and "Mismatch repair" were the most significantly enriched terms (Table S15). In the QCK 12 h vs. XCK 12 h comparison, 18,610 DEGs were assigned to 135 pathways in the KEGG database. "Plant-pathogen interaction", "Plant hormone signal transduction", "Monobactam biosynthesis", and "Circadian rhythm-plant" To deeply understand the DEGs between the two cultivars in response to waterlogging stress and reoxygenation, the DEGs were further investigated by performing a GO functional analysis. In the Q 0 h vs. X 0 h comparison, 6921 DEGs were assigned a GO term. Nine GO classes significantly fell into the category 'cellular component', 8 GO classes significantly fell into the category "molecular function", and 8 GO classes significantly fell into the category 'biological process' (Table S10). In the QCK 12 h vs. XCK 12 h comparison, 16,575 DEGs were assigned a GO term. Twenty-one GO classes significantly fell into the category "cellular component", 46 GO classes significantly fell into the category "molecular function", and 33 GO classes significantly fell into the category "biological process" (Table S11). In the QW 12 h vs. XW 12 h comparison, 7041 DEGs were assigned a GO term. Three GO classes significantly fell into the category "molecular function", and 13 GO classes significantly fell into the category "biological process" (Table S12). In the QCK 12 h + 2 h vs. XCK 12 h + 2 h comparison, 7374 DEGs were assigned a GO term. Fifteen GO classes significantly fell into the category "cellular component", 52 GO classes significantly fell into the category "molecular function", and 34 GO classes significantly fell into the category 'biological process' (Table S13). In the QW 12 h + 2 h vs. XW 12 h + 2 h comparison, 9443 DEGs were assigned a GO term. Seven GO classes significantly fell into the category "cellular component", 30 GO classes significantly fell into the category "molecular function", and 10 GO classes significantly fell into the category "biological process" (Table S14).
A KEGG pathway analysis of the abovementioned DEGs was also performed. In the Q 0 h vs. X 0 h comparison, 7623 DEGs were assigned to 135 pathways in the KEGG database. "DNA replication", "Isoquinoline alkaloid biosynthesis", "Tyrosine metabolism", and "Mismatch repair" were the most significantly enriched terms (Table S15). In the QCK 12 h vs. XCK 12 h comparison, 18,610 DEGs were assigned to 135 pathways in the KEGG database. "Plant-pathogen interaction", "Plant hormone signal transduction", "Monobactam biosynthesis", and "Circadian rhythm-plant" were the most significantly enriched terms (Table S16). In the QW 12 h vs. XW 12 h comparison, 7694 DEGs were assigned to 134 pathways in the KEGG database. "Mismatch repair", "Homologous recombination", "DNA replication" and "Nucleotide excision repair" were the most significantly enriched terms (Table S17). In the QCK 12 h + 2 h vs. XCK 12 h + 2 h comparison, 8000 DEGs were assigned to 135 pathways in the KEGG database. The primary terms included "Plant hormone signal transduction", "Plant-pathogen interaction", "Flavonoid biosynthesis" and "Isoflavonoid biosynthesis" (Table S18). In the QW 12 h + 2 h vs. XW 12 h + 2 h comparison, 10,028 DEGs were assigned to 135 pathways in the KEGG database. The primary terms included "Plant-pathogen interaction", "Regulation of autophagy", "Plant hormone signal transduction", and "Ubiquitin mediated proteolysis" (Table S19).

Hormone Response and Biosynthesis
Ethylene is a key hormone meditating the response to waterlogging in plants. In this study, two DEGs that were related to ethylene signaling and four ethylene synthesis genes were detected ( Figure 4). ACO (CL4948.Contig1), ACS1 (CL4253.Contig1), ERF1 (Unigene42398), ACS6 (CL15412.Contig4), and ACS7 (Unigene24656) were significantly up-regulated at 12 h of waterlogging and down-regulated after 2 h of recovery ( Figure 4). In addition, these genes were more induced in "Nannongxuefeng" than in "Qinglu", indicating that "Nannongxuefeng" had stronger ethylene synthesis activities, which is consistent with the ethylene measurement results (Table 1). These three ACSs were down-regulated after 2 h of recovery, but ethylene production increased, which may be due to the reoxygenation. ERF2 (CL11043.Contig1) was suppressed by waterlogging for 12 h and was up-regulated by reoxygenation for 2 h in "Nannongxuefeng", while ERF2 remained stable in the three samples of "Qinglu" (Figure 4), suggesting that ERF2 may play a different role in the two cultivars in waterlogging and recovery.

Other DEGs Involved in Waterlogging and Recovery in the Two Cultivars
Group VII ERFs are main transcription factors that regulate the expression of hypoxia-related genes and the response to waterlogging stress via the N-end rule pathway. Here, we analyzed five group VII ERFs and other N-end rule pathway genes ( Figure 5). Three group VII ERFs, HRE2 Contig4), and PCO2 (CL16190.Contig1) were up-regulated at 12 h of waterlogging stress and down-regulated after 2 h of reoxygenation in "Qinglu" and "Nannongxuefeng" (Figure 5). Except for PCO1, the other five gene transcripts were more abundant in "Nannongxuefeng" after waterlogging ( Figure 5), thus revealing stronger N-end rule pathway responses in "Nannongxuefeng" than in "Qinglu". However, another two VII ERFs, i.e., RAP2.12-like (Unigene52252 and CL7783.Contig3), were continuously down-regulated after waterlogging and recovery ( Figure 5). The different expression modes of the group VII ERFs indicated that these five members may perform different functions and have different mechanisms in the response to waterlogging. ROS signaling and changes in homeostasis are among the most vital responses to waterlogging. Here, we found eight DEGs that are related to the ROS signaling pathway ( Figure 5). Three PODs (CL15514.Contig3, Unigene15714, and Unigene10463) and AOX1a (CL9414.Contig3) were more up-regulated by waterlogging of 12 h in "Nannongxuefeng" than in "Qinglu" (Figure 5). Another three DEGs, i.e., CAT (CL3646.Contig2), CSD1 (CL8215.Contig1), and APX1 (CL5448.Contig2), and ROS signaling and changes in homeostasis are among the most vital responses to waterlogging. Here, we found eight DEGs that are related to the ROS signaling pathway ( Figure 5). Three PODs (CL15514.Contig3, Unigene15714, and Unigene10463) and AOX1a (CL9414.Contig3) were more up-regulated by waterlogging of 12 h in "Nannongxuefeng" than in "Qinglu" (Figure 5). Another three DEGs, i.e., CAT (CL3646.Contig2), CSD1 (CL8215.Contig1), and APX1 (CL5448.Contig2), and were more up-regulated during reoxygenation stage in "Nannongxuefeng" than in "Qinglu" ( Figure 5). Thus, the reactive oxygen scavenging enzymes in "Nannongxuefeng" had stronger activities than those in "Qinglu". In addition, the ROS generating gene RBohC (Unigene18044) was more up-regulated in "Qinglu" than in "Nannongxuefeng" during waterlogging stage ( Figure 5). The higher up-regulation of the ROS generating genes and lower responses of reactive oxygen scavenging enzymes may exacerbate the sensitivity to waterlogging in "Qinglu".
A heatmap analysis was also performed to study the DEGs that are related to anaerobic respiration and carbohydrate metabolism in the two different chrysanthemum cultivars during waterlogging ( Figure 5). Six key anaerobic respiration genes were found. Two PDCs (Unigene34986 and Unigene45721), three ADHs (Unigene15445, CL2541.Contig5, and Unigene6115), and SUS1 (CL20525.Contig1) were significantly up-regulated at 12 h of waterlogging and were down-regulated after 2 h of recovery in the two cultivars ( Figure 5). However, the three ADHs were more up-regulated in "Qinglu" than in "Nannongxuefeng" during the waterlogging stage ( Figure 5). The higher expression of the ADHs may accelerate the energy consumption during waterlogging, leading to sensitivity to waterlogging in "Qinglu".

Differential Transcription of Other Novel Genes Under the Waterlogging and Recovery Conditions
In addition to abovementioned DEGs, certain novel genes were identified in chrysanthemum under the waterlogging condition. Eight DEGs were further investigated ( Figure 6). Per1-like protein (CL18277.Contig2) and Snakin-2-like isoform X1 (Unigene878) were significantly up-regulated in "Nannongxuefeng" at 12 h of waterlogging, but down-regulated in "Qinglu" (Figure 6). Heat shock protein 83-like (Unigene31781), Chaperone protein ClpB1 (Unigene16908 and Unigene24844), autophagy-related protein 18a (ATG18a; Unigene420), and multiprotein-bridging factor 1c (MBF1c; Unigene15845) were more up-regulated in "Nannongxuefeng" than in "Qinglu" at 12 h of waterlogging ( Figure 6). In addition, Defensin SD2 (Unigene31384) was significantly up-regulated in "Nannongxuefeng" during reoxygenation ( Figure 6). However, the functions of these DEGs under the waterlogging and recovery conditions must be further investigated. were more up-regulated during reoxygenation stage in "Nannongxuefeng" than in "Qinglu" ( Figure  5). Thus, the reactive oxygen scavenging enzymes in "Nannongxuefeng" had stronger activities than those in "Qinglu". In addition, the ROS generating gene RBohC (Unigene18044) was more up-regulated in "Qinglu" than in "Nannongxuefeng" during waterlogging stage ( Figure 5). The higher up-regulation of the ROS generating genes and lower responses of reactive oxygen scavenging enzymes may exacerbate the sensitivity to waterlogging in "Qinglu". A heatmap analysis was also performed to study the DEGs that are related to anaerobic respiration and carbohydrate metabolism in the two different chrysanthemum cultivars during waterlogging ( Figure 5). Six key anaerobic respiration genes were found. Two PDCs (Unigene34986 and Unigene45721), three ADHs (Unigene15445, CL2541.Contig5, and Unigene6115), and SUS1 (CL20525.Contig1) were significantly up-regulated at 12 h of waterlogging and were down-regulated after 2 h of recovery in the two cultivars ( Figure 5). However, the three ADHs were more up-regulated in "Qinglu" than in "Nannongxuefeng" during the waterlogging stage ( Figure  5). The higher expression of the ADHs may accelerate the energy consumption during waterlogging, leading to sensitivity to waterlogging in "Qinglu".

Discussion
Water is an important environmental factor that is affecting plant growth, morphology, physiological and biochemical metabolism and geographical distribution. However, frequent waterlogging restrains the growth of plants and decreases yield, causing major disasters in agriculture [33]. Therefore, the identification of the plant waterlogging tolerance mechanism has always been an important worldwide agricultural problem. In the present study, we used RNA-seq technology to reveal the transcriptome changes in two chrysanthemum cultivars in response to waterlogging and post-stress reoxygenation. First, we investigated the morphological changes in "Nannongxuefeng" and "Qinglu" after waterlogging stress (Figure 1). After six days of waterlogging stress treatment, "Qinglu" showed a sensitive phenotype with poor growth potential, while "Nannongxuefeng" showed a waterlogging-tolerant phenotype and continued to grow well ( Figure 1). Thus, "Nannongxuefeng" is more tolerant than "Qinglu" to waterlogging.

Ethylene Results in Phenotypical Differences Between "Nannongxuefeng" and "Qinglu" after Waterlogging
Ethylene is a vital hormone that is induced by waterlogging stress [19,50]. In this study, the previously reported increase in ethylene production in chrysanthemum after waterlogging was confirmed [46]. We created four treatment groups and detected the ethylene changes during the waterlogging and reoxygenation stages. The waterlogging stress promoted ethylene production in "Nannongxuefeng" and "Qinglu" (Table 1). In addition, the ethylene synthesis precursor ACC combined with the waterlogging treatment promoted greater ethylene production than the waterlogging treatment alone, and the ethylene biosynthesis inhibitor AVG combined with the

Discussion
Water is an important environmental factor that is affecting plant growth, morphology, physiological and biochemical metabolism and geographical distribution. However, frequent waterlogging restrains the growth of plants and decreases yield, causing major disasters in agriculture [33]. Therefore, the identification of the plant waterlogging tolerance mechanism has always been an important worldwide agricultural problem. In the present study, we used RNA-seq technology to reveal the transcriptome changes in two chrysanthemum cultivars in response to waterlogging and post-stress reoxygenation. First, we investigated the morphological changes in "Nannongxuefeng" and "Qinglu" after waterlogging stress (Figure 1). After six days of waterlogging stress treatment, "Qinglu" showed a sensitive phenotype with poor growth potential, while "Nannongxuefeng" showed a waterlogging-tolerant phenotype and continued to grow well ( Figure 1). Thus, "Nannongxuefeng" is more tolerant than "Qinglu" to waterlogging.

Ethylene Results in Phenotypical Differences Between "Nannongxuefeng" and "Qinglu" after Waterlogging
Ethylene is a vital hormone that is induced by waterlogging stress [19,50]. In this study, the previously reported increase in ethylene production in chrysanthemum after waterlogging was confirmed [46]. We created four treatment groups and detected the ethylene changes during the waterlogging and reoxygenation stages. The waterlogging stress promoted ethylene production in "Nannongxuefeng" and "Qinglu" (Table 1). In addition, the ethylene synthesis precursor ACC combined with the waterlogging treatment promoted greater ethylene production than the waterlogging treatment alone, and the ethylene biosynthesis inhibitor AVG combined with the waterlogging treatment led to a lower ethylene production than the waterlogging treatment alone (Table 1). Interestingly, elevated ethylene production was also observed in both cultivars after two hours of reoxygenation, except for under the control treatment condition ( Table 1). The increase in the ethylene production after reoxygenation has been investigated in plants [49]. We speculate that the higher oxygen density results in an increase in ethylene production after reoxygenation. The ethylene production was higher in "Nannongxuefeng" than that in "Qinglu" ( Table 1), regardless of the treatment group. The striking ethylene synthesis response may be a vital adaption for coping with waterlogging and reoxygenation in the waterlogging-tolerant cultivar "Nannongxuefeng".

Changes in the Expression of Transcription Factors under the Waterlogging and Reoxygenation Conditions
Transcription factors are key regulators of the response to various biotic or abiotic stresses [51]. Here, 22 differentially expressed TFs from the AP2/ERF, bHLH, and MYB families were identified (Figure 3). Many studies have shown that these TF families are involved in abiotic stress and they positively improve plant tolerance [24,[52][53][54]. The N-end rule pathway plays a crucial role in the sensing of the oxygen density in plants [36,37,55]. The key DEGs in the N-end rule pathway were analyzed in this study. RAP2.12, RAP2.3, HRE2, ATE, PCO1, and PCO2 were up-regulated by waterlogging and down-regulated by reoxygenation in "Qinglu" and "Nannongxuefeng" (Figure 5). The key hypoxia regulators, i.e., RAP2.12, RAP2.3, and HRE2, in "Nannongxuefeng" were more highly induced than those in "Qinglu" (Figure 5), while in Arabidopsis, RAP2.12 and RAP2.3 are constitutively expressed [23]. In addition, two other group VII ERFs, i.e., RAP2.12-like (Unigene52252 and CL7783.Contig3), were down-regulated under the waterlogging condition ( Figure 5), suggesting that the functions performed by the group VII ERFs, which are the key regulators in the sensing of oxygen in the N-end rule pathway, may explain certain differences between chrysanthemum and Arabidopsis thaliana under hypoxic conditions.
Among the other TFs, MaRAP2-4 is a positive regulator of waterlogging tolerance via Mentha by regulating the energy-consuming processes and binding the AtSWEEET10 promoter in Arabidopsis [56]. TaMyb1 (Triticum aestivum Myb transcription factor 1) is strongly induced in the roots under the hypoxia condition in wheat [57]. HYPOXIA RESPONSE ATTENUATOR1 (HRA1) is a trihelix family TF that is highly up-regulated under hypoxic conditions and it acts as a negative regulator of the key hypoxia transcription factor RAP2.12 [58]. In addition, HRA1 negatively regulates the activation of its own promoter [58]. This regulatory mechanism maintains the balance of hypoxia-responsive genes, thereby preventing the excessive accumulation of specific proteins [58]. In this study, DREB2A, ERF017, ERF109 (Unigene30139 and Unigene22576), DREB1D, AIL5, and two MYB-related proteins were more strongly suppressed during the reoxygenation stage in "Nannongxuefeng" than "Qinglu" (Figure 3), suggesting that these TFs are negative regulators under reoxygenation conditions. ERF109 (Unigene27829), ERF (Unigene18230), bHLH93 (Unigene49514) , bHLH130 (Unigene35505 and CL2224.Contig1), SPATULA, MYB59, MYB86, and MYB5 were more strongly induced by waterlogging in "Nannongxuefeng" than in "Qinglu" (Figure 3). These differentially expressed TFs may be regulators that meditate the waterlogging response and stronger tolerance in "Nannongxuefeng". The functions of these TFs must be further investigated.

Waterlogging and Recovery Lead to Changes in Hormonal Responses and Biosynthesis-Related DEGs
The differential expression of genes that are related to ethylene response and signaling during waterlogging has been reported [59,60]. For example, in a comparative analysis of the transcriptome profiles of two cucumber genotypes ("Zaoer-N", which is waterlogging-tolerant, and "Pepino", which is waterlogging-sensitive) under waterlogging conditions, ACS and ACO were up-regulated in both lines, but the levels were higher in "Zaoer-N" [61]. Here, we selected pivotal genes in the ethylene synthesis pathway, i.e., ACS1, ACS6, ACS7, and ACO, and two pivotal genes for ethylene signaling, i.e., ERF1 and ERF2, for further investigation (Figure 4). ACS1, ACS6, ACS7, and ACO were more highly induced under waterlogging in "Nannongxuefeng" than in "Qinglu", suggesting that "Nannongxuefeng" has higher ethylene synthesis activities than "Qinglu" (Figure 4), which is consistent with the ethylene measurement results ( Table 1). The ACS genes were down-regulated during reoxygenation, while ethylene production was higher after reoxygenation. We speculate that these contradictory phenomena are due to the higher oxygen density caused by the reoxygenation, which promoted the transformation of the remaining ACC to ethylene [49].
Ethylene and IAA are key hormones that are meditating the root growth under waterlogging stress in plants. For example, the accumulation of ethylene promotes the formation of adventitious roots (AR) in waterlogged Rumex plants [21]. Higher accumulation of ethylene in flooded Rumex palustris Sm. increases the sensitivity of the root-forming tissues to endogenous IAA, thereby initiating the formation of AR [62]. The treatment with the ethylene biosynthesis inhibitor AVG and the auxin transport inhibitor 1-naphthylphthalamic acid (NPA) can result in a reduction of AR formation in waterlogged Solanum lycopersicum [63]. In this study, except for ethylene-related DEGs under waterlogging, several DEGs that are involved in the IAA response were also identified. ARF7 (Unigene43854), ARF19-like (Unigene41775) were more induced in "Qinglu" than in "Nannongxuefeng" during the waterlogging stage ( Figure 4). Auxin transporter-like protein 2 (Unigene17953) were more in "Nannongxuefeng" under waterlogging stress ( Figure 4). However, ARF7-like (Unigene38340) and ARF5 (Unigene31182) were more suppressed by waterlogging in "Nannongxuefeng" (Figure 4), suggesting that these DEGs may play different roles in Chrysanthemum roots under waterlogging stress. These ethylene and IAA-related DEGs may be the regulators that meditate the root development and AR formation during waterlogging. The functions of these DEGs must be further investigated.
In Arabidopsis, AtERF1 and AtERF2 are induced by various abiotic stresses, such as drought, salt, and cold [64]. AtERF1 and AtERF2 function as transcription activators and bind the GCC-box, thereby meditating down-stream gene expression [64]. In chrysanthemum, CmERF1 was induced by the waterlogging stress and was suppressed by reoxygenation ( Figure 4). However, waterlogging suppressed the expression of CmERF2 in the two chrysanthemum cultivars. In addition, CmERF2 was strongly induced during reoxygenation in "Nannongxuefeng" (Figure 4). Thus, CmERF1 and CmERF2 likely play different roles than AtERF1 and AtERF2 in chrysanthemum under waterlogging and reoxygenation conditions.

ROS Signaling Pathway, Anaerobic Respiration and Carbohydrate Metabolism
The up-regulation of ROS-scavenging enzyme encoding genes that are caused by waterlogging stress has been reported [70,71]. Here, three PODs and AOX1a were more significantly induced by waterlogging in "Nannongxuefeng" than in "Qinglu" (Figure 5), suggesting that "Nannongxuefeng" has higher antioxidant enzyme activities than "Qinglu" under waterlogging stress conditions. The level of RBohC for ROS generation was much higher under waterlogging in "Qinglu" than that in "Nannongxuefeng" (Figure 5), indicating the presence of higher ROS accumulation in "Qinglu". The higher ROS-scavenging enzyme activities and the lower ROS accumulation likely contribute to tolerance to waterlogging in "Nannongxuefeng".
Most environment-induced DNA damages originate in the disruption of ROS homeostasis [72], resulting in DNA single-and double-strand breaks [73][74][75]. For example, a ROS-scavenging enzymes double mutant apx1/ cat2, activates the DNA damage response and suffers more than wild type in Arabidopsis [76]. In this study, 'Mismatch repair', 'Homologous recombination', and 'Nucleotide excision repair' were significantly enriched in the KEGG analysis of the QW 12 h vs. XW 12 h comparison (Table S17). This result may be involved in the difference of ROS signaling responses between two cultivars. The higher ROS-scavenging enzyme activities during waterlogging may help to reduce the DNA damage in "Nannongxuefeng".
When hypoxic environments shift to normoxic conditions, significant changes can cause severe damages to plants [77], and ROS accumulation could also occur due to the rapid reoxygenation [13,78,79]. In this study, CAT, Cu/Zn, SOD, and APX were more significantly induced by reoxygenation in "Nannongxuefeng" than in "Qinglu", suggesting that certain ROS-scavenging enzymes play a role in the response to reoxygenation. Higher ROS-scavenging enzyme activities during reoxygenation may help chrysanthemum to cope with the damage that is caused by reoxygenation and adapt to environmental changes.
In plants, under hypoxic or waterlogging conditions, anaerobic respiration becomes the major form of ATP generation to cope with the energy crisis that is caused by the low oxygen level [5,7]. Transcriptome analyses of Cerasus sachalinensis roots in response to short-term waterlogging have been performed and indicated that the changes in the fermentation mechanism are the first coping mechanisms in response to hypoxic conditions [80]. Representative marker genes related to sucrose synthase and the glycolytic and fermentation pathways, such as PDC, ADH, and SUS, are up-regulated during the waterlogging stage [80]. In this study, PDC1, ADH, and SUS1 were significantly up-regulated by waterlogging and down-regulated by reoxygenation in the two chrysanthemum cultivars ( Figure 5). We speculate that the changes in the oxygen density may regulate the expression of these genes. Similar patterns in the expression of PDC1, ADH, and SUS1 were observed in "Nannongxuefeng" and "Qinglu", suggesting that the mechanism of anaerobic respiration and carbohydrate metabolism is highly conserved.

Other Genes Involved in Waterlogging and Reoxygenation
Several novel genes with differential expression patterns between the two cultivars under the waterlogging or reoxygenation conditions have also been investigated (Figure 6), including Defensin SD2, Snakin-2-like, HSP83-like, ClpB1, ATG18a, and MBF1c. Except for Defensin SD2, the other seven DEGs were more highly induced in "Nannongxuefeng" than in "Qinglu". The heat-inducible transcription factor HsfA2 is strongly induced by anoxia, and the over-expression of HsfA2 in Arabidopsis seedlings could enhance the tolerance to low oxygen [81]. Here, we found that the heat-related genes HSP83-like and MBF1c were more strongly induced by waterlogging in "Nannongxuefeng", but the function of these two genes and the crosstalk between heat and waterlogging stress in chrysanthemum must be further investigated.
In conclusion, we used RNA-seq technology to compare the global transcriptomes of two chrysanthemum cultivars in response to waterlogging and post-waterlogging reoxygenation. The DEGs between "Nannongxuefeng" and "Qinglu" varied greatly at the transcription level. In addition, "Nannongxuefeng" produces more ethylene than "Qinglu" under the waterlogging and reoxygenation conditions. Based on these results, a hypothetical model of the response to waterlogging and reoxygenation in chrysanthemum is presented (Figure 8). This study provides a comprehensive transcriptomic analysis for a better understanding of the response mechanism against waterlogging and reoxygenation, and more effective guidance for breeding chrysanthemum.

Plant Materials
The "Qinglu" and "Nannongxuefeng" Chrysanthemum cultivars were obtained from the Chrysanthemum Germplasm Resource Preserving Center of Nanjing Agricultural University, China. All of the plants were propagated with cuttings. The medium contained a 2:2:1 (v/v) mixture of perlite, vermiculite and leaf mold. The rooted seedlings were grown in a greenhouse under constant conditions of 22 °C during the day and a minimum of 15 °C during the night. The relative humidity varied from 70 to 75% (m/m), and no artificial light was provided. The healthy, uniform cuttings were transplanted into a small pot (4 cm down-diameter, 5 cm up-diameter, 9 cm deep) containing a 3:1:1 mixture of garden soil, perlite and vermiculite until the plants developed [10][11][12] leaves. Then, all of the plants were transferred to the culture room under the 16 h photoperiod (80-100 μmol/m 2 /s illumination) at 22 ± 1 °C conditions for waterlogging treatment.

Waterlogging Treatments
The waterlogging treatment was performed, as previously described by Yin et al. [46]. The pots of chrysanthemums were placed in a 28 cm × 14 cm × 14 cm container filled with tap water, and the water surface was 3 cm above the soil surface. The flooding depth was maintained on a daily basis. Control plants remained well-watered (60% soil moisture) throughout the experiment. The two cultivars were subjected to the waterlogging treatment for 0 d and 6 d for phenotype observation, and 0 h, 12 h, and 12 h + 2 h (2 h reoxygenation) for RNA-seq. Each treatment included at least three replications.

Ethylene Measurements
Four groups of treated chrysanthemum were assessed: control, waterlogging, waterlogging + ACC, and waterlogging + AVG. For waterlogging + ACC and waterlogging + AVG treatments, the pots were subjected in 200 μm ACC solution or 10 μm AVG solution. The treated "Qinglu" and "Nannongxuefeng" root samples were collected at 0 h, 3 h, 6 h, 12 h, and 24 h after waterlogging and after 2 h of reoxygenation after 12 h of waterlogging, and immediately transferred to 22-mL vials and maintained for 6 h at 25 °C. At the indicated time points, 1 mL of headspace air was removed using a syringe and was injected into a GC9790 Plus gas chromatograph (Fuli, Taizhou, China), and the ethylene content was measured, as previously described [82]. At least three independent duplicates were performed. The statistical significance of ethylene measurements data was

Plant Materials
The "Qinglu" and "Nannongxuefeng" Chrysanthemum cultivars were obtained from the Chrysanthemum Germplasm Resource Preserving Center of Nanjing Agricultural University, China. All of the plants were propagated with cuttings. The medium contained a 2:2:1 (v/v) mixture of perlite, vermiculite and leaf mold. The rooted seedlings were grown in a greenhouse under constant conditions of 22 • C during the day and a minimum of 15 • C during the night. The relative humidity varied from 70 to 75% (m/m), and no artificial light was provided. The healthy, uniform cuttings were transplanted into a small pot (4 cm down-diameter, 5 cm up-diameter, 9 cm deep) containing a 3:1:1 mixture of garden soil, perlite and vermiculite until the plants developed 10-12 leaves. Then, all of the plants were transferred to the culture room under the 16 h photoperiod (80-100 µmol/m 2 /s illumination) at 22 ± 1 • C conditions for waterlogging treatment.

Waterlogging Treatments
The waterlogging treatment was performed, as previously described by Yin et al. [46]. The pots of chrysanthemums were placed in a 28 cm × 14 cm × 14 cm container filled with tap water, and the water surface was 3 cm above the soil surface. The flooding depth was maintained on a daily basis. Control plants remained well-watered (60% soil moisture) throughout the experiment. The two cultivars were subjected to the waterlogging treatment for 0 d and 6 d for phenotype observation, and 0 h, 12 h, and 12 h + 2 h (2 h reoxygenation) for RNA-seq. Each treatment included at least three replications.

Ethylene Measurements
Four groups of treated chrysanthemum were assessed: control, waterlogging, waterlogging + ACC, and waterlogging + AVG. For waterlogging + ACC and waterlogging + AVG treatments, the pots were subjected in 200 µm ACC solution or 10 µm AVG solution. The treated "Qinglu" and "Nannongxuefeng" root samples were collected at 0 h, 3 h, 6 h, 12 h, and 24 h after waterlogging and after 2 h of reoxygenation after 12 h of waterlogging, and immediately transferred to 22-mL vials and maintained for 6 h at 25 • C. At the indicated time points, 1 mL of headspace air was removed using a syringe and was injected into a GC9790 Plus gas chromatograph (Fuli, Taizhou, China), and the ethylene content was measured, as previously described [82]. At least three independent duplicates were performed. The statistical significance of ethylene measurements data was calculated by One-Way ANOVA method using SPSS 23.0 (SPSS Inc., Chicago, IL, USA). Duncan's test was applied to assess the differences between treatments.

RNA Extraction, cDNA Library Construction and Sequencing
The abluent "Qinglu" and "Nannongxuefeng" root samples were immediately collected at 0 h, 12 h, and 12 h + 2 h, snap-frozen in liquid nitrogen and stored at −80 • C. Each sample was 0.5 g and three biological replicates were performed. The total RNA was isolated from each sample using the RNAiso Plus reagent (TaKaRa Bio, Tokyo, Japan) following the manufacturer's protocol. The quality and concentration of each RNA sample were verified using an ND-1000 Spectrophotometer (NanoDrop, Wilmington, DE, USA), and only RNA samples delivering the OD260/280 ratio of 1.8-2.1 and the OD260/230 ratio > 1.8 were retained. A 3 µg pool of RNA was formed by combining 1 µg from each biological replicate that was used for the construction of libraries.
The total RNA was treated with DNase I (TaKaRa Bio, Tokyo, Japan), and Oligos (dT) were used to isolate the mRNA. The mRNA was fragmented using a fragmentation buffer. Then, the first cDNA strand was synthesized using the mRNA fragments as templates, which was achieved by random hexamer priming. The second cDNA strand was synthesized by a reaction that was driven by DNA polymerase I (TaKaRa Bio, Tokyo, Japan). The short double-stranded DNA (dsDNA) fragments were purified using a QiaQuick PCR Extraction Kit (Qiagen, Hilden, Germany) and resolved with EB buffer for the end repair and single nucleotide A (adenine) addition. Sequencing adapters were ligated onto dsDNAs and suitable fragments were selected for the PCR amplification.

Transcriptome Data Processing and de novo Assembly
Clean reads were obtained from the raw reads after filtering the low-quality reads, adaptor-polluted reads, and reads with a high content of unknown bases (N). After the filtering, the remaining clean reads were stored in the FASTQ format. The de novo assembly using the clean reads was performed using Trinity to obtain the unigenes. Using BLAST+ (v2.2.23) software (http:// blast.ncbi.nlm.nih.gov/Blast.cgi) by default parameters, the unigenes were assigned functions based on the homologs present in the NR (www.ncbi.nlm.nih.gov/refseq), NT (www.ncbi.nlm.nih.gov/nuccore), Swiss-Prot (www.uniprot.org), KEGG (www.genome.jp/kegg/kegg1.html), COG (www.ncbi.nlm. nih.gov/COG/), and GO (geneontology.org) databases [83]. Using BLAST2GO (v2.5.0) software (https://www.blast2go.com) by default parameters with NR annotation to get the GO annotation [84]. The expression levels of the unigenes were calculated using the FPKM (fragments per kilobase of transcript per million fragments mapped) method [85]. The DEGs were defined according to an FDR ≤ 0.001 and a |log 2 fold change| ≥ 1, based on the PossionDis method [86]. The datasets of transcriptome are available in the NCBI repository, Accession No. for library SRP131844.

qRT-PCR Validation and Analysis
Eight DEGs were selected from all DEGs to validate the reliability of the libraries. The primers that were used for these validations were designed using Primer 5.0 software, and qRT-PCR was performed, as previously described by Ren et al. [87]. The chrysanthemum EF1α gene (GenBank accession number KF305681) was used as a reference, and the gene primers are listed in Table S20. Three biological replicates were performed per sample.