The S. tonkinensis seeds are derived from Pingxiang, Jiangxi Province, China, and plug seedlings were cultivated in a greenhouse at the Xiashu Forest Farm, Zhenjiang, Jiangsu Province, China (118°40′E, 31°41′N). The young plug seedlings with an average height of 15 cm were transplanted in permeable non-woven containers containing a medium of organic matter: perlite: vermiculite (v/v: 6:2:2) to continuous growth through common water and nutrient management. The waterlogging treatments were started at 8 a.m. on 5 September 2021 when seedlings reached the average height of 55 cm. The seedlings were planted plot-in-plot (Supplementary Figure 1A). The inside basin is the permeable non-woven container, and the outside basin is a water-proof lotus plot. During the experimental period, the daily average temperature varied between 25.9 and 29.6°C. Relative humidity (RH) ranged from 62 to 76%. The daylight varied between 12.6 and 12.7 h with 13,316 KJ/m² average daily solar irradiation. The waterlogging treatment lasted for 4 days, and the water surface elevated the matrix by 3–5 cm. A total of 16 seedlings were destructively collected for each sample counting the survival rate every day. The primary root samples (Supplementary Figure 1B) were collected at 0, 9, 24, and 72 h after waterlogging treatment, immediately frozen in liquid nitrogen and kept at –80°C refrigerators for further processing. The sample without waterlogging treatment (0 h) was set as the control group. A number of three biological replicates for each analytical experiment of the transcriptome analyses and the quantitative real-time PCR (qRT-PCR) experiments were harvested.

The content of three phytohormones including salicylic acid (SA), abscisic acid (ABA), and indole-3-acetic acid (IAA) was determined. The endogenous phytohormones were extracted using acetonitrile solution. About 0.5 g of fresh root samples was first mixed with acetonitrile solution at a volume ratio of 1:10, followed by centrifugation, purification, and nitrogen blowing. The measurements were taken using the AGILENT 1,290 high-performance liquid chromatography (HPLC) (AGILENT Company, Palo Alto, CA, United States) series SCIEX-6500Qtrap mass spectrometer (MS/MS) (AB Company, Framingham, MA, United States). Chromatographic separation of five types of phytohormones was performed using a Poroshell 120 SB-C18 reverse-phase chromatography column (2.1 mm × 150 mm, 2.7 μm). Phase A (methanol/0.1% formic acid solution) and phase B (water/0.1% formic acid solution) of HPLC-MS/MS were utilized to conduct gradient elution at a constant flow rate of 0.3 ml/min at 30°C. In this experiment, multiple reaction detection (MRM) scan mode was applied to detect the substance. For the substances to be tested, 2 or more fragment ions were determined. The ratio of their peak time and response value corresponded to the standard substance was analyzed, to identify the detected objects. The standard substances that include D-ABA, D-IAA, and D-SA were purchased from Sigma-Aldrich (Shanghai, China) to draw the standard curve to qualify the detected objects.

Approximately 0.2 g of fresh root samples was collected to determine 1-aminocyclopropane-1-carboxylic acid (ACC). Under hypoxia conditions, plants can generate the ET precursor (ACC) in roots and vertically transfer it to the xylem. Subsequently, ACC would be oxidized to ET in shoots. ACC would be determined to replace ET due to that ET is volatile and hard to be detected (Banga et al., 2010; Bradford et al., 2018). It was isolated by 4°C precooling deionized water for 2 h. After centrifugation and purification, 2 μl of samples was displaced into AGILENT 1,290 high-performance liquid chromatography series 6420A mass spectrometer (AGILENT Company, Palo Alto, CA, United States) for further detection. Poroshell 120 SB-C18 reverse-phase chromatography column (2.1 mm × 150 mm, 2.7 μm) was also applied. Phase A (acetonitrile solution) and phase B (water/0.1% formic acid solution) of HPLC-MS/MS were utilized to conduct gradient elution at 35°C. The ACC standard was bought from J&K Scientific Ltd. (Beijing, China).

Root samples collected at 0, 9, 24, and 72 h were used for transcriptome analyses. Total RNAs were extracted using the Plant RNA Kit (Ambion) following the manufacturer’s protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). The samples with RNA integrity number (RIN) ≥ 7 were subjected to the subsequent analysis. The libraries were constructed using TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, United States) according to the manufacturer’s instructions. Then, these libraries were sequenced on the Illumina Hiseq 4000 Sequencing platform, and 125/150-bp paired-end reads were generated.

The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China) using Illumina Hiseq 4000 Sequencing platform. Raw data (raw reads) were processed using Trimmomatic (Bolger et al., 2014). The reads containing poly-N and the low-quality reads were removed to obtain the clean reads. After removing adaptor and low-quality sequences, the clean reads were assembled into expressed sequence tag clusters (contigs) and de novo assembled into the transcript using Trinity (Grabherr et al., 2013) in the paired-end method. The longest transcript was chosen as a unigene based on the similarity and length of a sequence for subsequent analysis.

The function of the unigenes was annotated by alignment of the unigenes with the NCBI non-redundant (NR), SwissProt, and clusters of orthologous groups for eukaryotic complete genome (KOG) databases using BLASTx (Altschul et al., 1990) with a threshold E-value of 10^–5. The proteins with the highest hits to the unigenes were used to assign functional annotations thereto. Based on the SwissProt annotation, Gene Ontology (GO) classification was performed by the mapping relation between SwissProt and GO term. The unigenes were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) database to annotate their potential metabolic pathways.

The DESeq2 method was applied to standardize the count number of each sample gene, and the base mean value was used to estimate the expression level, which was calculated as fragments per kilobase per million mapped reads (FPKM). Besides, the multiple of difference was calculated, and negative binomial (NB) distribution test was used to test the significance of the difference. Finally, the differential protein-coding genes were screened according to the different multiple and different significance test results (Love et al., 2014). Expression levels of unigenes were calculated. The DEGs between different time points were identified with p < 0.05 and | log2FC| > 1 (Anders and Huber, 2010).

To validate the RNA-seq results, a total of eight transcripts ([IAA], PYL, EBF1_2, EIN3, MYC2, NPR1, AHP, and AUX1/LAX), which were associated with phytohormones biosynthesis at each of sampling points, were verified by qRT-PCR. The primers for each DEG are displayed in Table 1. All reactions were done on a StepOne Real-Time PCR System using SYBR Green Dye (Applied Biosystems, Foster City, CA, United States; Takara, Dalian, China). Relative gene expression was evaluated using the 2^–ΔΔCt method with 18S ribosomal RNA as an internal control.

The data analysis included a basic descriptive analysis followed by an analysis of variance (ANOVA). Duncan and Pearson R correlation tests were performed using SPSS 23.0 for Windows (SPSS Science, Chicago, IL, United States). The p-values less than 0.05 and 0.01 were both considered to indicate significance between groups. For the elaboration of graphs, Origin 2018 (OriginLab, Northampton, MA, United States) was used.

1-Aminocyclopropane-1-carboxylic acid and three phytohormones at 0, 9, 24, and 72 h were quantified by the HPLC-MS/MS method. The group of 0 h was regarded as the control group to compare with waterlogged groups. For each time point, it had three biological replicates. ACC, a precursor of ET, was sharply raised to 369.2 ± 38.35 Aa ng/g⋅FW after 72 h of waterlogging stress. During the first day of waterlogging, a slight variation in ACC content in roots appeared among three sampling points (Figure 1A). Compared to the control group (0 h), the ABA content in roots under waterlogging treatment decreased, especially for the seedlings that experienced 9 and 24 h of waterlogging. Nonetheless, ABA content was pointing to a bit of a rebound on the third day (12.88 ± 0.75 Bb ng/g⋅FW), which had significant differences from other groups (Figure 1B). The IAA content was elevated substantially after 9 h first and dropped subsequently on the third day (Figure 1C). Waterlogging inhibited SA accumulation in S. tonkinensis roots. Compared to the seedlings without waterlogging, SA content descended by approximately 30% during the first day of waterlogging and rapidly dropped by 76% with the aggravation of flooding stress on the third day (Figure 1D).

A total of 78.73 G of clean data were obtained after transcriptome sequencing of 12 samples was completed. The effective data amount of each sample ranged from 6.12 to 6.92 G, Q30 bases ranged from 94.65 to 95.51%, and the average GC content was 46.12% (Table 2). A total of 84,601 unigenes were spliced with a total length of 81,389,823 bp and an average length of 962.04 bp (Table 3). Concerning the mapping statistics, the total reads of each sample ranged from 43,206,426 to 48,601,426 with a ratio of 85.53–89.67% when comparing reads to unigenes (Table 4). The 12 samples from 0, 9, 24, and 72 h were analyzed through the DESeq2 method to determine FPKM (Figure 2). For the 0-h group, the gene number of FPKM ranged from 1 to 10 was 33,537, 32,958, and 36,850 for each sample, respectively. 0h3 sample differenced with 0h1 and 0h2 samples. In the 9-h group, the gene number of FPKM over 10 and FPKM ranged from 1 to 10 was also similar between each sample.

The quality control was applied after data output statistics via box-whisker plot and the purpose of principal component analysis (PCA). In the transcriptome results of S. tonkinensis seedlings under waterlogging stress, phytohormone-related unigenes were selected to conduct further analysis. The related unigenes were renewedly managed to draw Supplementary Figure 2. It was indicated the distribution of gene expression levels in Supplementary Figure 2A. The median gene expression amount for each group and the height of these boxes were basically at the same level. Even though these samples originated from the same tissue (primary roots) of S. tonkinensis seedlings, three biological replicates in each group still had a slight gene expression difference. PCA is to assess the biological repeatability of samples within groups and differences between groups. As shown in Supplementary Figure 2B, the long distance in the graph between each group illustrated obvious group differences. It was remarkable that the 24h1 sample was close to the 72-h group, manifesting that the 24h1 sample was more susceptible to waterlogging stress. For each group, the biological repeatability performed decently as the distance in the graph between each sample in a group was close.

In this section, total DEGs and DEGs related to plant hormones in 9 vs. 0 h, 24 vs. 0 h, and 72 vs. 0 h were analyzed via GO classification, respectively (Supplementary Table 2). Massive GO pathways were filtered to the top 30 of enrichment for the following analysis. In the 9 vs. 0-h comparison, a total of 4,637 DEGs were assigned to three main GO functional categories, including biological process (3,752, 80.91%), cellular component (4,046, 87.25%), and molecular function (4,032, 86.95%). The biological process, cellular component, and molecular function were independently assigned into 10 sub-categories with the majority of unigenes in “cellular process” (up 68.24%, down 66.21%), “organelle” (up 58.84%, down 61.75%), and “binding” (up 59.52%, down 54.20%). As shown in Figures 3B,C and a total of 167 DEGs related to plant hormones were separately annotated to biological process (155, 92.81%), cellular component (125, 74.85%), and molecular function (143, 85.63%). In the biological process, the upregulated unigenes in “auxin-activated signaling pathway” and “ethylene biosynthetic process” ranked 2nd and 5th with 14 and 6 unigenes, respectively. A total of 12 downregulated unigenes in “auxin-activated signaling pathway” ranked top 1.

In 24 vs. 0-h and 72 vs. 0-h comparison, a total of 8,238 DEGs and 7,146 DEGs were, respectively, assigned into three main GO functional categories. The most abundant sub-categories of biological process, cellular component, and molecular function 24 vs. 0 h and 72 vs. 0 h were the same as that in 9 vs. 0 h (Figures 4A, 5A). Concerning the biological process, “response to stimulus” with 2,865 upregulated unigenes (34.78%) and 2,131 unigenes (25.87%) were displayed in the top 10 biological processes in 24 vs. 0 h (Figure 4A). In 72 vs. 0 h, the smaller number of the upregulated unigenes (2,228, 31.18%) involved in “response to stimulus” is displayed in Figure 5A when compared to that in 24 vs. 0 h. In plant hormone signal transduction and biosynthesis pathway, 253 EDGs in 24 vs. 0 h and 227 DEGs in 72 vs. 0 h were assigned into three main GO functional categories (Figures 4, 5). Generally, the number of downregulated DEGs related to the “auxin-activated signaling pathway” was increased to 22 unigenes with the elongation of waterlogging. Additionally, the upregulated DEGs in the “ethylene-activated signaling pathway” played an essential role in the biological process during the first day of waterlogging. The upregulated DEGs involved in the “abscisic acid-activated signaling pathway” remained high until the third day of waterlogging.

Through KEGG pathway analysis, the unigenes related to the biosynthesis of ET, ABA, IAA, and SA were acquired in Figure 6. S-adenosylmethionine synthetase (SAMS, [EC:2.5.1.6]), 1-aminocyclopropane-1-carboxylate synthase (ACS, [EC:4.4.1.14]), and aminocyclopropanecarboxylate oxidase (ACO, [EC:1.14.17.4]) are the three core enzymes in the biosynthesis of ET in plants. We identified 47 unigenes encoding SAMS, 4 encoding ACS, and 7 unigenes encoding ACO. Except for several SAMS genes, the expression patterns of these unigenes mostly peaked in 9 h after waterlogging treatment. A number of three key enzymes including 9-cis-epoxycarotenoid dioxygenase (NCED, [EC:1.13.11.51]), xanthoxin dehydrogenase (ABA2, [EC:1.1.1.288]), and abscisic-aldehyde oxidase (AAO3, [EC:1.2.3.14]) are actively participated in ABA biosynthesis, which were encoded by 12 unigenes, 1 unigene and 2 unigenes, respectively. FPKM of NCED raised from 0.42–0.96 to 1.94–3.96 after experiencing 3-day waterlogging. The expression patterns of these unigenes showed a down–up trend. A number of two core enzymes, which are L-tryptophan—pyruvate aminotransferase (TAA1, [EC:2.6.1.99]) and indole-3-pyruvate monooxygenase (YUCCA, [EC:1.14.13.168]), play an essential role prior to IAA biosynthesis. FPKM of TAA1 dropped from 1.61 ± 0.33a to 0.14 ± 0.08b after 3-day waterlogging. S. tonkinensis seedlings carried out phenylalanine ammonia-lyase (PAL, [EC:4.3.1.24]) to generate SA in roots. We determined 6 unigenes encoding PAL with a downward trend after 3 days under waterlogging stress.

After DEG acquisition, KEGG pathway significance analysis was conducted to classify these DEGs (Kanehisa et al., 2008). In the pathway of phytohormone signal transduction, a total of 48 positively regulated unigenes were identified and were targeted by 18 phytohormone-associated genes (ABF, 3 ARR-A, 2 BAK1, BIN2, 2 CTR1, EIN3, 3 ERF1, 4 ETR/ERS, 4 GH3, 10 [IAA], JAZ, PP2C, PR1, PYL, 5 SAUR, 5 SnRK2, 2 TCH4, and TGA) when comparing 9-h group and 0-h group. However, 33 unigenes were downregulated and were classified into 17 types of unigenes (ABF, 2 AHP, 2 ARR-B, 4 AUX1/LAX, BSK, DELLA, ETR/ERS, 5 [IAA], JAZ, 2 K14486/ARF, MYC2, 2 PP2C, 2 PR1, PYL, SAUR, 3 SnRK2, and 3 TGA) (Supplementary Figure 3 and Supplementary Table 1). A total of 11 DEGs and 14 DEGs in the auxin signaling transduction pathway were upregulated and downregulated, respectively (Supplementary Figure 6A). Compared to 0 h, 55 DEGs were upregulated and 25 genes were determined (2 ABF, 2 ARR-A, BAK1, BSK, BIN2, 2 CTR1, CYCD3, EIN3, 3 ERF1, 3 ETR/ERS, 3 GH3, GID1, 8 [IAA], JAR1_4_6, 2 JAZ, K14486/ARF, MYC2, 2 PP2C, 3 PR1, 2 PYL, 7 SAUR, 5 SnRK2, TCH4, TGA, and TIR1) in roots which had been waterlogged for 24 h, whereas 50 negatively regulated unigenes were sorted into 23 classifications associated with plant hormones (ABF, 4 AHP, 2 ARR-A, 2 AUX1/LAX, BRI1, BSK, 3 CYCD3, DELLA, EBF1_2, ETR/ERS, 2 GH3, 7 [IAA], 2 JAZ, 3 K14486/ARF, MYC2, NPR1, 4 PP2C, 2 PR1, 2 PYL, 3 SAUR, 3 SnRK2, TCH4, and TGA) (Supplementary Figure 4). Compared to the control group (0 h), the 72-h group owned 32 upregulated DEGs and 63 downregulated DEGs (Supplementary Figure 4). These DEGs were identified and sorted to 17 positively related classifications (2 ABF, BAK1, 2 BSK, EIN3, 2 ERF1, ETR/ERS, 3 GH3, GID1, [IAA], JAZ, K14486/ARF, 2 PP2C, 3 PR1, 2 PYL, 7 SAUR, SnRK2, and TGA) and 27 negatively related classifications (ABF, AHK2_3_4, 6 AHP, 4 ARR-A, 2 ARR-B, 4 AUX1/LAX, BAK1, BSK, COI-1, 3 CYCD3, DELLA, 2 EBF1_2, ERF1, ETR/ERS, 2 GH3, 9 [IAA], 2 JAZ, 4 K14486/ARF, MYC2, NPR1, 2 PP2C, 2 PR1, PYL, 6 SAUR, 2 SnRK2, TCH4, and TGA) (Supplementary Figure 5 and Supplementary Table 1). Generally, it was indicated that the number of significantly downregulated unigenes was increased with the time extension of waterlogging, especially in the auxin signal transduction pathway (Supplementary Figure 6 and Supplementary Table 1).

Based on the results of our transcriptomic analysis, eight genes known to be involved in phytohormones signal transduction were selected to confirm their mRNA expression levels by qRT-PCR. The results of qRT-PCR are displayed in Supplementary Figure 7, which were consistent with their variation trends in transcriptional profiles.

Ethylene is an essential phytohormone, which performs as a chemical messenger in many physiological processes throughout a plant life cycle (Schaller and Voesenek, 2015). It also gets involved in a wide range of combating responses to abiotic and biotic stresses (Müller and Munné-Bosch, 2015). Oxygen deprivation around roots caused by waterlogging induces ET accumulation in most plant species (Voesenek et al., 2003; Jackson, 2008; Alpuerto et al., 2016). Under waterlogging conditions, ET improves plant adaption to flooding by triggering the formation of adventitious roots in tomatoes (Vidoz et al., 2016), playing an important role in the ROS-signaling pathway and antioxidant system (Khan et al., 2020). Coudert et al. (2010) demonstrated that ET facilitated the formation of adventitious roots in rice. In the process of ET biosynthesis, methionine is catalyzed by SAMS to synthesize s-adenosylmethionine (SAM). Subsequently, SAM is carboxylated by ACS into ACC, followed by the oxidization of ACC to form ET under ACO regulation (Sauter et al., 2013; Chen et al., 2016). As indicated in Figure 6, the expression level of ACO (E1.14.17.4) peaked at 9 h after waterlogging, followed by a downward trend with the limited O₂. It was suspected that the ACC oxidation process was enforced at the early stage to generate ET, whereas ACC started to massively accumulate under hypoxia conditions with the low expression of ACO (Figure 1A).

Commonly, the ET receptor (ETR/ERS) can stimulate the kinase activity of constitutive triple-response 1 (CTR1) without ET. The activated CTR1 will repress the activity of ET insensitive 3 (EIN3) transcription factors (Guo and Ecker, 2003; Alonso and Stepanova, 2004; Chen et al., 2005). In our study, both CTR1 and ETR/ERS kept the same trend of obvious upregulation throughout waterlogging process (Supplementary Table 1). Previous studies identified two F box proteins which were called EBF1 and EBF2 in Arabidopsis, respectively. EBF1 and EBF2 mediate ubiquitination or proteasome pathway to negatively regulate ET signal transduction. EBF1_2 also called EIN3-binding F-box protein leads to degradation of EIN3 when EBF1_2 is overexpressed, consequently causing plants insensitive to ET (Guo and Ecker, 2003; Potuschak et al., 2003), which is the opposite to the results in Supplementary Table 1. EIN3 is an important positive regulator of the downstream ET pathway. The stabilization of EIN3 also refers to the ability to rapidly sense hypoxia signaling and activate the ET pathway in the early stage (Chao et al., 1997; Xie et al., 2015). This discovery was consistent with the results that DEGs related to the “ethylene-activated signaling pathway” were induced under waterlogging during the first day in Figures 3, 4. Genetic analysis proved that ET response factor 1 (ERF1) belongs to the AP2/ERF superfamily and participates downstream of EIN3 and other ET pathway components. ERF1, a positive ET-responsive gene, also can be regulated by EIN3 (Chao et al., 1997; Ma et al., 2014). In our study, EIN3 was positively regulated since the downregulated EBF1_2 cannot suppress the expression of EIN3, further activating the role of ERF1 in the downstream “ethylene-activated signaling pathway” (Guo and Ecker, 2003; Alonso and Stepanova, 2004; Chen et al., 2005). In general, the upregulated EIN3 indicated the ability to rapidly sense hypoxia signaling and activate the ET pathway in S. tonkinensis seedlings at the early stage of waterlogging (Chao et al., 1997; Xie et al., 2015).

Abscisic acid is a critical plant hormone for plant growth, which can accumulate in plants in response to adversity by improving turgor pressure (Vishwakarma et al., 2017). However, the results revealed that ABA accumulation in the roots of S. tonkinensis seedlings has been inhibited (Figure 1B). During ABA biosynthetic procedure, 9-cis-violaxanthin is prepared to produce xanthoxin under the catalysis of NCED. Then, ABA2 induces the dehydrogenation of xanthoxin to form abscisic aldehyde, which is eventually oxidized to synthesize ABA under AAO3 catalysis (Nambara and Poll, 2005). The gene expression of ABA biosynthesis-related unigenes took on a down–up trend (Figure 6), which was consistent with the dynamic pattern of ABA content in Figure 1B.

Abscisic acid can be perceived by its receptor proteins from the pyrabactin resistance 1-like (PYL) family, further sequestrating type 2C protein phosphatases (PP2C) while activating sucrose non-fermenting 1-related protein 2 (SnRK2). Subsequently, ABA signaling involves phosphorylation of the downstream ABRE-binding factor (ABF), which recognizes and regulates the expression of ABA-responsive genes, thereby enhancing the resistance of plants at expense of growth (Fujii et al., 2009; Ma et al., 2009; Zhao Y. et al., 2018; Belda-Palazón et al., 2020; Cui et al., 2020). As shown in Supplementary Table 1, PYL was negatively regulated, leading to weak suppression of PP2C and negatively regulated SnRK2 in roots of S. tonkinensis seedlings. It was inferred that upregulated ABF could recognize and regulate ABA in roots of S. tonkinensis seedlings under flooding stress (Mustilli et al., 2002; Yoshida et al., 2006; Fujii et al., 2007; Umezawa et al., 2009; Vlad et al., 2009).

Auxins (e.g., IAA) are known to be the major phytohormones involved in the process of the formation of adventitious roots, which is a vital adaptive trait under waterlogging stress (Vidoz et al., 2010; Zhang et al., 2015). Morris (1993) demonstrated that IAA played an important role in regulating plant growth and development to cope with hypoxia. The variation of IAA in roots of S. tonkinensis had an upward trend first and a downward trend subsequently on the third day (Figure 1C). The IAA biosynthesis process is a part of tryptophan metabolism. TAA1 uses tryptophan to form indole pyruvic acid (Stepanova et al., 2008). Mashiguchi et al. (2011) reported that the YUC gene family catalyzes the direct conversion of indole-3-pyruvic acid monohydrate to IAA. Won et al. (2011) discovered that TAA1 and YUC genes play a synergistic role in auxin biosynthesis. The expression level of unigenes related to IAA biosynthesis presented the fluctuation. Generally, FPKM of TAA1 and YUCCA at 72 h after waterlogging was decreased when compared to the control group, indicating the declined content of IAA in the 72-h group (Figure 1C).

Previous studies indicated that auxin-responsive genes (Aux/IAA, GH3, and SAUR) and auxin response factor genes (ARFs) participated in the process of the “auxin-activated signaling pathway” to strengthen plant resistance to abiotic stress (Hagen and Guilfoyle, 2002; Xie et al., 2020). ARFs are combined with Aux/IAA transcriptional repressors when auxin concentration is low, thereby inhibiting auxin-responsive genes (GH3 and SAUR). With the increase of auxin, transport inhibitor response 1 (TIR1) starts to recruit Aux/IAA, contributing to the ubiquitination and degradation of Aux/IAA (Dharmasiri et al., 2005; Mockaitis and Estelle, 2008; Maraschin et al., 2009). Subsequently, ARFs are released with proteolysis, activating auxin transcription (Leyser, 2018). Li et al. (2020) also reported that SAUR72, ARF6, and [IAA] genes were significantly upregulated whereas the GH3 gene was significantly downregulated in waterlogged plants. Nonetheless, Zhang et al. (2016) discovered that several auxin-related genes were negatively regulated in leaves of waterlogged cotton. Further, Kramer and Ackelsberg (2015) found that ARF17 was positively correlated with GH3 in the leaves and roots of Populus deltoides “DHY,” which was opposite to the consequences in Supplementary Table 1. In the first two groups, SAUR, GH3, and [IAA] genes were upregulated whereas K14486/ARF and Auxin influx carrier (AUX1/LAX) were downregulated. AUX1/LAX plays an essential role in the polar auxin transport (PAT), which is responsible for transporting auxin from root tips to the elongation zone (Péret et al., 2012).

Salicylic acid, a chemical messenger, can be synthesized in abundance in plants to improve plant tolerance under abiotic and biotic stresses. It is a β-hydroxy phenolic acid, which induces growth, development, and biochemical traits of plants at even minuscule concentrations (Arif et al., 2020). Miura and Tada (2014) expounded that SA physiologically got involved in plants by inducing stomatal closure and regulating plant morphology. The shikimic acid pathway is one of the synthetic routes for SA. Phenylalanine originated from shikimic acid and is catalyzed by PAL to synthesize trans-cinnamic acid (trans-CA) (Mishra and Baek, 2021). Subsequently, trans-CA is converted to o-coumaric acid (o-CA), and o-CA is directly converted to SA (León et al., 1995). The expression trend of PAL was identical to the dynamic pattern of SA content (Figure 1D).

In this study, three SA-responsive genes (PR1, TGA, and NPR1) participated in SA signaling pathway. Johnson et al. (2003) revealed that transcription factor TGA (TGA) stimulated the pathogenesis-related protein 1 (PR1) promoter in vivo in the SA and regulatory protein NPR1 (NPR1)-dependent manner. Besides, the gene expression of PR is regulated by NPR1 with transcription factors such as TGA (Zhou et al., 2000; Kim and Delaney, 2002). Ali et al. (2018) claimed that PR genes get significantly altered under adversities. This was consistent with our results that PR1 was significantly altered between waterlogging samples and the control group (0 h) (Supplementary Table 1).

The plant perception of abiotic stress triggers the crosstalk between phytohormones in response to adversities. This crosstalk forms a signaling network (Harrison, 2012). In the study, ACC showed an upward trend whereas ABA, IAA, and SA took on a downward trend from 0 to 72 h under waterlogging stress. Yamauchi et al. (2020) disclosed that AUX/IAA-mediated auxin signaling possibly plays a role in ET-dependent aerenchyma formation in rice roots under hypoxia conditions. As an alternative to the accumulation of ABA, ET may act as a root-sourced messenger transported as ACC in water stress circumstances (Jackson, 1997). Sharp (2002) proposed that endogenous ABA inhibited ET generation under anaerobic conditions. The interaction between ABA and ET contributes to the maintenance of shoot and root growth (Sharp, 2002). Yin et al. (2015) pointed out that ET and ABA seem to act synergistically or antagonistically to control plant growth and development. Benschop et al. (2005) and Saika et al. (2007) expounded that accumulated ET was not only an obstacle to NCED expression but also caused a breakdown of ABA into phaseic acid, thereby decreasing ABA. Table 5 indicates that waterlogging had an opposite effect on the content of ET precursor (ACC) and ABA (r = −0.973). Previous research found ABA accumulation to foster stomatal closure (Savchenko et al., 2014). On the contrary, Chini et al. (2007) and Thines et al. (2007) explained that auxins regulated stomatal opening positively. ABA can cooperate with auxin to facilitate plant growth, owing to that ABA treatment improves the gene expression of ARF2 but decreases IAA7 to activate the auxin signaling pathway (Song et al., 2009; Wang et al., 2011). In this study, K14486/ARF was downregulated, but [IAA] was upregulated with the reduction of ABA content (Figure 1B and Supplementary Table 1). Combined with the decreased IAA content, it was concluded that waterlogging had a similar impact on the variation of ABA and IAA content (r = 0.982) in S. tonkinensis (Table 5).

Overall, the upregulated EIN3 and downregulated EBF1_2 assisted in ethylene signaling at an early stage and kept S. tonkinensis seedlings sensitive to ET. Nonetheless, ET biosynthesis was inhibited by the low expression level of ACO on the third day of waterlogging. In addition, the “ET activated signaling pathway” failed to form adventitious roots and aerenchyma in S. tonkinensis seedlings (Figure 7). The “auxin-activated signaling pathway” was enriched in both the upregulated and downregulated genes in 9 vs. 0 h and 24 vs. 0 h, whereas it was only enriched in the downregulated genes in 72 vs. 0 h. The falling expression level of PAL resulted in the low accumulation of SA content at 72 h (Figure 7). Moreover, the decreased transcription of PYL weakened the ABA signal transduction pathway in the roots of S. tonkinensis seedlings under waterlogging stress.

One limitation in our study is the large biological variability since these seedlings originated from seeds. In the future, spraying exogenous ABA and SA on the seedlings from cutting propagation will be conducted to reinforce the suspected conclusion. Also, further exploration of the interaction of these genes plays an essential role in the understanding of the crosstalk between phytohormones and theoretically supports the molecular mechanism of S. tonkinensis seedlings’ resistance to waterlogging.