Xiaoxiao Zhang
Xiang Liu3†
Yonghong Hu- 1College of Landscape Architecture and Arts, Northwest A&F University, Yangling, Shaanxi, China
- 2Shanghai Key Laboratory of Plant Functional Genomics and Resources, Shanghai Chenshan Botanical Garden, Shanghai, China
- 3School of Ecological Technology and Engineering, Shanghai Institute of Technology, Shanghai, China
Paeonia ostii, a widely cultivated tree peony species in China, is a resourceful plant with medicinal, ornamental and oil value. However, fleshy roots lead to a low tolerance to waterlogging in P. ostii. In this study, P. ostii roots were sequenced using a hybrid approach combining single-molecule real-time and next-generation sequencing platforms to understand the molecular mechanism underlying the response to this sequentially waterlogging stress, the normal growth, waterlogging treatment (WT), and waterlogging recovery treatment (WRT). Our results indicated that the strategy of P. ostii, in response to WT, was a hypoxic resting syndrome, wherein the glycolysis and fermentation processes were accelerated to maintain energy levels and the tricarboxylic acid cycle was inhibited. P. ostii enhanced waterlogging tolerance by reducing the uptake of nitrate and water from the soil. Moreover, transcription factors, such as AP2/EREBP, WRKY, MYB, and NAC, played essential roles in response to WT and WRT. They were all induced in response to the WT condition, while the decreasing expression levels were observed under the WRT condition. Our results contribute to understanding the defense mechanisms against waterlogging stress in P. ostii.
Introduction
Paeonia ostii, belonging to section Moutan DC. of the genus Paeonia L. (Paeoniaceae), is a perennial woody shrub native to China (Zhang et al., 2017a). Paeonia ostii has been cultivated as ornamental and medicinal in China for more than 1500 years (Li et al., 2011; Hu and Han, 2018). In recent years, P. ostii was planted extensively to harvest the roots for the traditional Chinese medicine “Mudanpi” with functions of clearing heat and cooling blood, while recently, farmers have also started collecting the seeds for oil production (Hu and Han, 2018). It has been found that the seeds of P. ostii are rich in unsaturated fatty acids, especially α-linolenic acid (Li et al., 2015; Yu et al., 2016), which are beneficial for human health (Lefevre et al., 2004; Shahidi and Miraliakbari, 2005). Seed oil has been authenticated as a new resource of functional food since 2011 (Zhang et al., 2017b). Moreover, the area under cultivation of P. ostii is the largest among all the ten tree peony species and nearly 1500 cultivars in China (Zhang et al., 2018a), with the flower being white and monopetalous. In the middle and lower reaches of the Yangtze River, an irregular spatial and temporal distribution of precipitation frequently occurs, which results in water levels that exceed the requirement for P. ostii, causing waterlogging in the field (Ding et al., 2020). Therefore, P. ostii often fails to grow well in this region, and the main reasons for this phenomenon are waterlogging and high temperature (Wu et al., 2009; Hu et al., 2017).
Waterlogging can be defined as the saturation of soils with water, and it often occurs during the rainy season. Waterlogging is an adverse abiotic stress that limits the oxygen diffusion into the soil and creates a low oxygen (hypoxia) environment around the plant roots (Juntawong et al., 2014). Hypoxia can accelerate the anaerobic respiration, reduce the root activity, and result in energy shortage. Plants can temporarily maintain energy production to some extent during hypoxia caused by waterlogging, via glycolysis and ethanol fermentation (Pan et al., 2021). The genes encoding enzymes involved in these pathways such as pyruvate decarboxylase (PDC), alcohol dehydrogenase (ADH), and lactate dehydrogenase (LDH) were consistently upregulated under waterlogging stress (Yin et al., 2019). Besides, a key feature for the acclimation to hypoxia is to activate genes encoding enzymes involved in transcription regulation and signaling pathways in order to allow biological and physiological adjustments to the hypoxia conditions (Bailey-Serres et al., 2012). The ethylene response factor represented the highest number of significantly expressed TFs under hypoxic conditions followed by bHLH, MYB, NAC, and WRKY in soybean (Chen et al., 2016). As waterlogging recedes, these plant roots, which had adjusted to the reduced light and oxygen in murky water, are suddenly re-exposed to aerial conditions. The shift to an intensely reoxygenated environment poses an additional stress for the plant, namely oxidative stress. Therefore, waterlogging acts as sequential stress manipulating plant growth (Yeung et al., 2018). Our previous study showed that the physiological response to waterlogging of P. ostii was rapid and sensitive. Despite its failure to tolerate long-term waterlogging, P. ostii possesses a strong recovery ability in the critical frame (within 72 h after waterlogging) (Hu et al., 2017). However, the molecular mechanisms underlying the responses to waterlogging and waterlogging recovery in P. ostii remain largely unexplored.
As a next-generation sequencing (NGS) technique, RNA-sequencing has become an indispensable tool for transcriptome-wide analysis of differential gene expression and gene regulatory networks (Hu et al., 2020). RNA-seq provides a precise and comprehensive analysis of RNA transcripts that affect gene expression (Chang et al., 2019). It has recently been widely used for researching P. ostii gene expression in response to various stresses, including chill-induced stress (Gai et al., 2013), copper stress (Wang et al., 2016), and drought-induced stress (Zhao et al., 2019). However, for the species without reference genomes, NGS was insufficient to analyze their gene expression because of inaccurate reconstruction and lack of reliable expression estimation of transcript variants due to the short, sequenced reads (Steijger et al., 2013). Single-molecule real-time (SMRT) sequencing, a third-generation technique developed by Pacific Biosciences (PacBio) to sequence cDNA, can produce longer reads than NGS and avoids further assembly (Sharon et al., 2013). The combination of SMRT and NGS provided researchers with an efficient and effective way to determine gene expression patterns, thereby exploring unanswered biological questions, particularly for species lacking genome-related information (Minoche et al., 2015; Deng et al., 2018; Jia et al., 2018).
Plant roots are the tissues in direct contact with water under waterlogging stress. Thus, we hypothesized that identifying the genes upregulated and downregulated under hypoxia in P. ostii roots can deepen our understanding of the molecular regulatory pathways in response to waterlogging and waterlogging recovery. In the present study, the seedings with 3 d of waterlogging treatment and followed by a 7-d recovery treatment, were selected as subjects, according to our previous results (Hu et al., 2017). We constructed a de novo full-length transcriptomic database for seedlings roots under all conditions using the SMRT sequencing. Furthermore, RNA-seq was adopted to identify the differentially expressed genes (DEGs) and analyze the specific pathways involved in waterlogging and waterlogging recovery treatment. Thus, the results of this study will deepen our understanding of the defense mechanisms against waterlogging stress in P. ostii and provide theoretical support for the molecular breeding of P. ostii to withstand waterlogging stress.
Materials and methods
Plant materials
In the present study, 3-year-old potted P. ostii were used as study materials. On March 1, the well-grown plants cultivated in the field were transplanted into the pots, and then the plants were transported to the phytotron. The day and night temperatures were set to (22 ± 1)°C and (16 ± 1)°C, respectively, with a photoperiod of 16 h light/8 h dark, 60% relative humidity, 40% soil water content, and light intensity of 600 μmol m-2 s-1. On May 1, uniform plants were selected and randomly divided into three groups. In the control (CK) group for waterlogging treatment (WT), the sampling time was 0 h post-waterlogging. The WT water level was set at 2 cm above the soil, and the sampling time was 3 d after waterlogging. For the waterlogging recovery treatment (WRT), the soil surface water was removed after 3 d of waterlogging, and the sampling time was 7 d after recovery. Then, the same plants from the WT group were also used as the control for the WRT group. For RNA-seq, two biological replicates were performed, and each biological sample consisted of the roots of three uniform plants. The mixed samples from these three sampling points were collected for full-length sequencing. All these samples were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. Total RNA was extracted from the roots using TIANGEN RNA Prep Pure Plant Plus kit (Tiangen Biotech Co. Ltd., Beijing, China), and the quality and quantity of RNA were assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).
Measurement of morphological and physiological indices
The roots of plants from three treatments were washed with deionized water, and root morphology was photographed by Epson Perfection V700 Photo (Epson (China) Co., Ltd., China). Root tips were cut by a knife blade, fixed in FAA, and then conducted by using a Saffron-O and Fast Green Stain Kit (Solarbio, Beijing, China) based on the manufacturers’ instructions. The cell morphology of root tip was viewed with an optical microscope (BX43, Olympus, Tokyo, Japan). Root activity was detected by triphenyl tetrazolium chloride (TTC) reduction method (Hu et al., 2017). Relative electrical conductivity (REC) of the root was measured as described previously (Zhao et al., 2019). All tests were performed with three biological replicates.
Illumina transcriptome library preparation and sequencing
A total of 5 μg RNA per sample was used as input material to generate sequencing libraries using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out in the NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. The library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA) to identify cDNA fragments of about 300 bp length. Then, 3 μl USER Enzyme (NEB, USA) was incubated with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C. Then PCR was performed with Phusion® High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified with AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq 2500 platform, and paired-end reads were generated.
PacBio SMRT bell library preparation and sequencing
The Iso-seq template was prepared based on the protocol of Iso-Seq Template Preparation for Sequel Systems. First strand cDNA was synthesized from 800–1000 ng total RNA using Clontech SMARTer™ PCR cDNA Synthesis Kit (Clontech Laboratories, Inc., USA). The CDS Primer IIA was first annealed to the polyA+ tail of transcripts, followed by first-strand synthesis with SMARTScribe™ Reverse Transcriptase. Then, a sufficient amount of double-stranded cDNA was produced using large-scale PCR with Clontech PrimeSTAR GXL DNA Polymerase and 5’PCR Primer IIA (5’- AAGCAGTGGTATCAACGCAGAGTAC-3’). After PCR amplification, the products were selected using the BluePippin Size Selection System with the following bins for each sample: 1–2, 2–3, and > 3 kb. The DNA was repaired by DNA Damage Repair Mix (PacBio) and End Repair Mix (PacBio). Blunt Adapter was ligated to the cDNA using the ligase from PacBio. Exonucleases Exo III and EXO VII were added to the ligated cDNA, to reduce the amount of unrepaired DNA or linear DNA without a blunt adapter. Finally, the obtained cDNAs were measured using Qubit HS (Life Technologies, USA) and Agilent Bioanalyzer 2100. Sequencing reactions were performed using the PacBio Sequel sequencer (BGI-Shenzhen, China) with Sequel Sequencing Kit 2.1 and Sequel SMRT Cell 1M v2 Tray.
Data processing of PacBio sequencing reads
Raw sequencing data (also called raw polymerase reads) produced by the Pacific Biosciences Sequel system were processed following the IsoSeq protocol through the SMRT analysis package version 2.3.0 (Pacific Biosciences, https://www.pacb.com/products-and-services/analytical-software/smrt-analysis). These data have been deposited in the Genome Sequence Archive (https://bigd.big.ac.cn/gsa) in the National Genomics Data Center, under the accession number: CRA004521. The raw polymerase reads that had full passes >0 and the predicted consensus accuracy >0.75 were selected for producing ROIs (Reads of Insert). ROIs with a minimum length of 300 bp were classified into full-length and non-full-length transcript sequences, based on whether the 5’ primer, 3’ primer, and poly-A tail were all observed. The full-length sequences were processed to de novo consensus isoforms using the Iterative Clustering for Error Correction algorithm and then polished via the Quvier quality-aware algorithm. The de novo consensus isoforms with high quality (the expected Quiver accuracy ≥ 0.95) from each library were merged, and redundancy was removed using Cluster Database at High Identity with Tolerance (CD-HIT) (Fu et al., 2012) based on the sequence similarity, to obtain final unique full-length isoforms. The coding sequences of the isoforms were identified by Transdecoder (v3.0.1), and then the longest one was selected perform BLAST (Basic Local Alignment Search Tool).
Functional annotation
Final full-length isoforms were mapped to the NCBI non-redundant protein sequences (NR), NCBI non-redundant nucleotide sequence (NT), SwissProt (a manually annotated and reviewed protein sequence database), Kyoto Encyclopedia of Genes and Genomes (version 59) (KEGG), and Clusters of Eukaryotic Orthologous Groups (KOG) database by Blast software (version 2.2.23, https://blast.ncbi.nlm.nih.gov/Blast.cgi) (Altschul et al., 1990) with default parameters (under a threshold E-value ≤10-5) to get the isoform annotations. Gene Ontology (GO) annotations and functional classifications were obtained using the Blast2GO program (version 2.5.0, E-value ≤10-5, https://www.blast2go.com) (Conesa et al., 2005) based on NR annotations. Hmmscan software (version 5.11- 51.0) (Jones et al., 2014) as used to obtain the annotations from the Pfam database.
Illumina data analysis
The sequencing data was filtered with SOAPnuke (v1.5.2, https://github.com/BGI-flexlab/SOAPnuke) (Li et al., 2008), and clean reads were obtained and stored in FASTQ format. The clean reads were mapped to reference full-length transcriptome using HISAT2 (v2.0.4, http://www.ccb.jhu.edu/software/hisat/index.shtml) (Kim et al., 2015). Bowtie2 (v2.2.5, http://bowtiebio.sourceforge.net/%20Bowtie2%20/index.shtml) (Langmead and Salzberg, 2012) was applied to align the clean reads to the reference coding gene set and then the expression levels of genes were calculated using RSEM (v1.2.12, https://github.com/deweylab/RSEM) (Li and Dewey, 2011). Essentially, differential expression analysis was performed using the DESeq2 (v1.4.5, http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) with Q value ≤ 0.05. DEGs were identified using DESeq2 (Love et al., 2014) with Q value (adjust P value) < 0.001 and fold change ≥ 2 or ≤ -2. The identified DEGs were subsequently subjected to GO and KEGG enrichment using Phyper in the R package, with default Q value <= 0.05.
Quantitative real-time PCR (qRT-PCR)
First-strand cDNA was prepared from 1 μg of total RNA per sample, using a FastKing RT Kit with gDNase (Tiangen Biotech Co. Ltd., Beijing, China). Specific primers were designed for each of the 14 selected isoforms (Table S1). PCRs were performed on an ABI StepOnePlus® Real-Time PCR System (Applied Biosystems, California, USA), following the manufacturer’s instructions. Each reaction mixture (20 μl) contained 10 μl of TB Green Premix Ex Taq II (Tli RNaseH Plus) (Takara), 0.8 μl of each primer (10 μM), 0.3 μl of cDNA template (1 μg), and 8.1 μl of RNase-free water. PCR for each gene was performed in triplicate, with the following thermal cycling conditions: 95°C for 30 s; 40 cycles of 95°C for 5 s and 64°C for 30s; and 95°C for 15 s. Primer specificity was confirmed via melt curve analysis. The relative expression levels of the tested genes were calculated via the 2-ΔΔCt method, using the β-actin gene as internal controls.
Statistical analysis
The results were expressed as the mean ± standard deviation (SD). One-way analysis of variance (ANOVA) and Duncan’s multiple range tests were used to analyze the significance at a p level of 0.05 by SPSS software (version 19.0 for Windows; SPSS Inc., 2010). The gene expression heatmap was performed using the OmicShare tools (https://www.omicshare.com/tools).
Results
Morphological and physiological indices
As shown in Figure 1A, the fibrous roots of plants in CK group were numerous and white. After waterlogging, almost half of the fibrous roots fell off, and the remaining fibrous roots became brown. The tips of some main roots were rotten. For the WRT group, although the rot of main root was serious, more fibrous roots had developed. These new fibrous roots were white, and the root system was basically the same as that of CK group. In terms of the anatomy of root tip cell (Figure 1B), the size became large after waterlogging, but the cytoplasm in most cells decreased. For the WRT group, the cell size of root tip decreased, and the cytoplasm in most cells increased.
Figure 1 Morphological changes of the control (CK), waterlogging treatment (WT), and waterlogging recovery treatment (WRT) of P. ostii. and anatomic (A), Root morphology; (B), Cell morphology of root tip.
The plant in the CK group had high root activity (Figure 2A). After waterlogging, the root activity decreased significantly, and the activity level was about one sixth of that of CK group. When the soil surface water was removed, the root activity slightly recovered. Although the root activity of plants in WRT group had reached about twice that of WT group, it was still significantly lower than that of CK group. An opposite trend was observed in the REC of plant root under different treatments (Figure 2B). There was significant variation in the REC of plant root among three treatment groups, and the highest value was found in the root of WT group.
Figure 2 Physiological changes of the control (CK), waterlogging treatment (WT), and waterlogging recovery treatment (WRT) of P. ostii. (A), Root activity; (B), Relative electrical conductivity of root. Values represent mean ± standard deviation (SD), and letters indicate significant differences according to Duncan’s multiple range test (p < 0.05).
PacBio ISO-seq analysis
PacBio ISO-seq was used in this study to sequence the mixed root samples and to construct a de novo full-length transcriptomic database due to the lack of a P. ostii reference genome. Moreover, 11 cDNA libraries, including insert fragments of 1–2 kb, 2–3 kb, and >3 kb, were prepared. In summary, 1,272,956 polymerase reads were generated by the PacBio Sequel platform (Table S2). Polymerase reads with full passes > 0 and the predicted consensus accuracy > 0.75 were selected for producing ROIs, and 1,055,458 ROIs were obtained in all libraries (Table S3). Furthermore, 546,928 full-length non-chimeric (FLNC) reads with poly (A) tail signals, 5’ adaptor sequences, and 3’ adapter sequences were obtained (Table S4). According to the isoform-level clustering algorithm, the full-length sequences were processed to de novo consensus isoforms. The isoforms were polished using Quvier quality-aware algorithm. The results from each library were merged, and redundancy was removed using CD-HIT. Finally, the transcripts contained 187,564 unique full-length isoforms, and those with an average length of 2323 bp were considered in the reference transcriptome and used for further analysis (Table S5). We further evaluated the completeness of the assemblies by comparing them against a set of conserved plant genes in the Benchmarking Universal Single-Copy Orthologs (BUSCO, embryophyta_odb9 dataset) using BUSCO v2.0 pipeline. The complete sequence accounted for 89.77% of all sequences (Figure S1A).
Annotation of the full-length reference transcriptome
All the isoforms of the reference transcriptome were aligned to the public databases including NCBI non-redundant proteins (NR), NCBI non-redundant nucleotides (NT), SwissProt, Kyoto Encyclopedia of Genes and Genomes (KEGG), KOG (Clusters of Eukaryotic Orthologous Groups), Pfam, and Gene Ontology (GO) (Table 1). Our results indicated that a total of 156,402 isoforms (83.39%) were annotated in at least one database, and 48,392 isoforms (25.80%) had significant matches in all seven databases. The best-harbored isoforms in the reference transcriptome were in the NR database (147,828, 78.81%), while only 85,754 isoforms (45.72%) got hits in the Pfam database. Overall, 59,824 isoforms got hits in five databases (NR, KEGG, KOG, SwissProt, and Pfam) (Figure S1B).
Table 1 Information of function annotation.
For GO analysis, the isoforms were assigned to 54 functional groups, which were allocated into three ontologies: biological process (167,903 isoforms), cellular component (230,396 isoforms), and molecular function (146,236 isoforms) (Figure S1C). For the biological process category, isoforms involved in cellular processes (GO:0009987, 50,400 isoforms) and metabolic processes (GO:0008152, 45,217 isoforms) were highly represented. Cell (GO:0005623, 44,564 isoforms), cell part (GO:0044464, 43,542 isoforms), membrane (GO:0016020, 36,727 isoforms), organelle (GO:0043226, 33,368 isoforms), and membrane part (GO:0044425, 33,251 isoforms) were the five most functional terms in the cellular component ontology. Within the molecular function category, binding (GO:0005488, 67,324 isoforms) was the largest group, followed by catalytic activity (GO:0003824, 63,247).
KEGG pathways database contains a systematic analysis of inner-cell metabolic pathways and functions of gene products. Metabolic pathway analysis of unique isoforms was also carried out using the KEGG annotation system. The results showed that 118,052 isoforms were classified into five main categories (Figure S1D). Most of isoforms included in the metabolism category were involved in many pathways, such as global and overview maps, and carbohydrate metabolism, amino acid metabolism, and energy metabolism. The second highly enriched category was genetic information processing, and the isoforms were mainly associated with translation and fold, sorting, and degradation. There were only 4976 isoforms included in environmental information processing and focused on the environmental adaptation pathways.
Moreover, all the isoforms were subjected to a search in the KOG database for functional prediction and classification. Overall, 119,842 isoforms were functionally classified into 25 specific categories (Figure S1E). Among them, general functional prediction only (35905 isoforms) was the largest category, followed by signal transduction mechanisms (17490 isoforms), post-translational modification, protein turnover, chaperone (12900 isoforms), function unknown (9480 isoforms), transcription (8785 isoforms), and carbohydrate transport and metabolism (8608 isoforms). Only 175 isoforms were assigned to cell motility.
RNA-seq analysis
RNA samples from roots were individually prepared from non-waterlogging, waterlogging, and waterlogging recovery group plants and then sequenced on an Illumina instrument, to obtain comprehensive waterlogging and waterlogging recovery responses in P. ostii. Illumina sequencing data have been deposited in the Genome Sequence Archive under the accession number: CRA004511. In total, 377,050,830 raw reads were generated. After removing adaptor sequences, ambiguous nucleotides, and low-quality sequences, 337,972,250 clean reads were recorded. The lowest Q20 percentage and GC percentages were 96.89% and 44.51%, respectively, for the six libraries. The filtered clean reads were mapped to the reference transcriptome generated by PacBio ISO-seq. Of all the reads, 57.04–65.30% were mapped to the reference transcriptome, and 10.80–15.06% were uniquely mapped reads (Table 2).
Table 2 Information of RNA-seq raw data and clean reads mapped with the reference transcriptome.
Identification and functional profiles of DEGs
The expression levels of all isoforms were calculated using RSEM software and shown as FPKM values. DEGs were identified using DESeq2 with Q value (adjust P value) < 0.001 and fold change ≥ 2 or ≤ -2. A total of 2951 DEGs were obtained in response to waterlogging treatment (CK VS WT); 1358 DEGs were significantly upregulated, and 1593 DEGs were significantly downregulated (Figure 3A). As for the waterlogging recovery treatment (WT VS WRT), 1730 DEGs were observed, with 759 of them upregulated and 971 downregulated (Figure 3B).
Figure 3 The distribution and expression levels of DEGs in P. ostii roots. (A), Volcano plot of waterlogging treatment (CK VS WT); (B), Volcano plot of waterlogging recovery treatment (WT VS WRT); (C), Venn diagram of DEGs under waterlogging and its recovery treatment. The x-axis represents the log2 [Fold Change] values under the mean normalized expression of all isoforms (y-axis).
The DEGs were then subjected to the Gene Ontology (GO) analysis to determine their functional classification. Under the WT, 2951 DEGs could be categorized into 45 functional groups on three main categories (biological process, cellular component, and molecular function) (Figure 4A). As for the WRT, these three main categories included 47 groups (Figure 4B). In both treatments, cellular process (GO:0009987), and metabolic process (GO:0008152) were the two largest groups enriched in the biological process; cell (GO:0005623), membrane (GO:0016020), membrane part (GO:0044425), and organelle (GO:0043226) were the four most prominent groups enriched in the cellular component; catalytic activity (GO:0003824) and binding (GO:0005488) were the two largest groups enriched in the molecular function. Using enrichment analysis, the top 20 significant enriched GO terms identified in the DEGs of two treatments (Q value ≤ 0.05) were shown in (Figure S2). Oxidoreductase activity (GO:0016491) was the most enriched term in both WT and WRT. In addition, there were five identical terms: heme-binding (GO:0020037), iron ion binding (GO:0005506), tetrapyrrole binding (GO:0046906), oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (GO:0016705), and cofactor binding (GO:0048037) in Top 20 enriched terms of two treatments.
Figure 4 GO enrichment classification of DEGs in P. ostii roots. (A), Waterlogging treatment (CK VS WT); (B), Waterlogging recovery treatment (WT VS WRT).
The DEGs were also subjected to KEGG analysis to investigate their biological functions further; 1404 DEGs of CK VS WT were assigned to 124 pathways, and 855 DEGs of WT VS WRT were categorized into 119 pathways. Enrichment analysis for DEGs showed that 15 pathways and 11 pathways were significantly enriched (Q value ≤ 0.05) under the WT (Figure 5A) and WRT (Figure 5B), respectively. Furthermore, 8 pathways, tyrosine metabolism (ko00350), fatty acid degradation (ko00071), alpha-linolenic acid metabolism (ko00592), glycolysis/gluconeogenesis (ko00010), monoterpenoid biosynthesis (ko00902), isoflavonoid biosynthesis (ko00943), biosynthesis of amino acids (ko01230), and carbon fixation in photosynthetic organisms (ko00710), were significantly enriched in both treatment groups. Interestingly, the proportion of upregulated and downregulated genes in these pathways behaved oppositely in the two treatment groups. For example, in the glycolysis/gluconeogenesis pathways, the number of upregulated genes was more than that of downregulated genes in the WT group. However, the number of downregulated genes was more than that of upregulated genes in the WRT group.
Figure 5 KEGG enrichment analysis of DEGs in P. ostii roots. (A), Waterlogging treatment (CK VS WT); (B), Waterlogging recovery treatment (WT VS WRT).
Since transcription factors (TFs) could mediate the expression of genes involved in the waterlogging and waterlogging recovery, we analyzed the DEGs coding TFs. A total of 112 TFs were classified into 22 families based on their assigned protein in the WT group, which accounted for 3.80% of the DEGs. Of these, 80 TFs were upregulated, and 32 TFs were downregulated (Figures 6A, C). 76 TFs were classified into 18 families in the WRT group, which accounted for 4.40% of the DEGs. Among them, 19 TFs were upregulated and 57 TFs were downregulated (Figures 6B, D). Apetala 2/ethylene-responsive element binding protein (AP2/EREBP) and MYB family represented the most significantly expressed TFs (22) in the WT group, while the TFs belonging to AP2/EREBP (16) were the most significantly expressed in the WRT group, followed by the MYB family (13). On the whole, most of TFs were upregulated in the WT group but downregulated in the WRT group.
Figure 6 DEGs coding transcription factors based on their assigned protein families in P. ostii roots. (A), Waterlogging treatment (CK VS WT); (B), Waterlogging recovery treatment (WT VS WRT); (C), Heatmap of DEGs coding transcription factors under waterlogging treatment; (D), Heatmap of DEGs coding transcription factors under waterlogging recovery treatment.
Analysis of genes involved in key pathways
The pathways of glycolysis, fermentation, citrate cycle, ethylene biosynthesis, and nitrogen metabolism were shown in Figure 7, and the RNA-Seq expression of transcripts encoding the enzymes involved in these pathways varied under WT and WRT. There were 839 DEGs in both the WT and the WRT (Figure 3C), and 14 DEGs that were mainly involved in glycolysis, fermentation, citrate cycle, ethylene biosynthesis, nitrogen metabolism, water absorption, and transcriptional regulation were selected for further qRT-PCR analysis (Figure 8). These DEGs encoded ADH, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), PDC, pyruvate kinase (PK), isocitrate dehydrogenase (IDH), malate dehydrogenase (MDH), glutamate decarboxylase (GAD), aquaporin, nitrate transporter (NRT), 1-aminocyclopropane-1-carboxylate oxidase (ACO), WRKY, MYB, NAC, and AP2/EREBP, respectively. The expression levels of all fourteen isoforms significantly varied between the waterlogging and waterlogging recovery groups. Our results showed that IDH, MDH, aquaporin, and NRT were downregulated in the WT group, but upregulated in the WRT group. However, the changes of the other ten isoforms were opposite in the WT and WRT groups. Expression analysis of all the selected isoforms of a few genes was carried out using qRT-PCR to validate the sequencing results, which showed similar patterns to those observed in the RNA-Seq data.
Figure 7 Waterlogging and waterlogging recovery caused genes encoding proteins involved in glycolysis, fermentation, citrate cycle (A), ethylene biosynthesis (B), nitrogen metabolism (C). The red character represents DEGs in both CK VS WT and WT VS WRT. INV, invertase; SUS, sucrose synthase; G6PC, glucose-6-phosphatase; PGM, phosphoglucomutase; GPI, glucose-6-phosphate isomerase; PFK, 6-phosphofructokinase; FBA, fructose-bisphosphate aldolase; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; PGAM, phosphoglycerate mutase; PPH, phosphopyruvate hydratase; PK, pyruvate kinase; PDC, pyruvate decarboxylase; ADH, alcohol dehydrogenase; LDH, lactate dehydrogenase; CS, citrate synthase; MDH, malate dehydrogenase; SDH, succinate dehydrogenase; SCS, succinyl-CoA synthetase; OGOR, 2-oxoglutarate/2-oxoacid ferredoxin oxidoreductase; IDH, isocitrate dehydrogenase; GDH, glutamate decarboxylase; ACS, 1-aminocyclopropane-1-carboxylate synthase; ACO, 1-aminocyclopropane-1-carboxylate oxidase; AP2/EREBP, apetala 2/ethylene-responsive element binding protein; NTR, Nitrate transporter; NR, Nitrate reductase; NiR, Nitrite reductase; NOR, Nitric-oxide reductase; NOS, Nitrous-oxide reductase.
Figure 8 Quantitative real-time PCR validation of key genes. Actin was used as the reference gene. Expression values were normalized such that the expression levels of CK were set to 1. GAPDH, glycer-aldehyde-3-phosphate dehydrogenase; PK, pyruvate kinase; PDC, pyruvate decarboxylase; ADH, alcohol dehydrogenase; MDH, malate dehydrogenase; IDH, isocitrate dehydrogenase; GDH, glutamate decarboxylase; ACO, 1-aminocyclopropane-1-carboxylate oxidase; AQP, aquaporin; NRT, nitrate transporter. Values represent mean ± standard deviation (SD), and letters indicate significant differences according to Duncan’s multiple range test (p < 0.05).
Discussion
Roots are the first organ to sense waterlogging, and hence play a key role in the waterlogging stress response (Langan et al., 2022). The root activity of P. ostii under waterlogging decreased significantly in comparison with the control, and increased as the waterlogging receded. Once a plant detects waterlogging, its priority is to reinstate the oxygen supply to the roots, which can be achieved by altering its root morphology and anatomy (Pedersen et al., 2021). These changes were also found in P. ostii. The tips of main roots were rotten, and a lot of fibrous roots fell off under waterlogging. Anatomic analysis showed that the cell size became larger due to a large amount of water absorption. This change in cell size may be due to the damage of cell membrane permeability. REC can reflect the degree of cell membrane damage (Zhao et al., 2019). After waterlogging, REC increased significantly, and then slightly decreased when the soil surface water was removed. These results indicated that waterlogging stress destroyed the cell membrane structure, and caused membrane permeability to become large, which was also supported by the anatomy data.
With advancements in technology, the studies on plant response to adverse abiotic stress are no longer limited to detecting the physiological characteristics but also identifying their molecular mechanisms. RNA-Seq has been successfully used to elucidate the response of plants to various environmental stresses, such as cold (Pang et al., 2013), salt (Postnikova et al., 2013), and drought (Zhao et al., 2019). Moreover, the combined hybrid approach of RNA-Seq and SMRT sequencing provides a more accurate integrated analysis of transcriptomes than using NGS alone (Chao et al., 2019). The clear advantage was that the traditional gene cloning to obtain full-length cDNA sequences, which is expensive, time-consuming, and inefficient, was not required (Deng et al., 2018). In the present study, we constructed the complete and accurate transcriptome of P. ostii roots using PacBio ISO-seq. Our study is novel because this is also the first full-length transcriptome sequencing of P. ostii.
In summary, we constructed three libraries (inserted size 1–2 kb, 2–3 kb, and > 3 kb) and obtained 4.71, 4.6, and 4.49 GB of raw data, respectively. A total of 13.8 GB of raw data was generated, and 187,564 full-length isoforms with an average length of 2323 bp were mapped to seven databases (NR, NT, KEGG, GO, KOG, SwissProt and Pfam). These data were similar to the full-length transcriptome profiling for flower bud of P. suffruticosa ‘High Noon’ (Chang et al., 2019) and we highlight how the gene sequence information would be instrumental for the further study of tree peony. The mRNAs of six P. ostii roots, belonging to CK, WT, and WRT groups, were sequenced using RNA-Seq and referred to as the transcriptome generated by ISO-Seq. We identified 2951 DEGs and 1730 DEGs in response to WT and WRT, respectively, and 839 DEGs were common between the two treatment groups. Our study demonstrated that the mechanisms involved in the WRT group were the reverse of the WT group and had some unique response mechanisms. Moreover, we identified response mechanisms involved in the WT group that did not manifest in the WRT group. Previous studies solely focused on waterlogging (Lee et al., 2014; Ren et al., 2017) or waterlogging recovery (Yeung et al., 2018). However, in our study, the gene expression differences underlying these two adverse stress conditions were investigated simultaneously, with the aim to deepen the comprehensive understanding of this sequential stress.
In the condition of waterlogging stress, plants have two response mechanisms. Hypoxia escape syndrome enables plants to overcome waterlogging through the elongation of petioles, enables stems, and leaves to get more oxygen through the formation of large lenticels, adventitious roots, and aerenchyma facilitates oxygen transportation and storage. However, these processes require much higher ATP and may lead to their death (McDonald et al., 2002; Glenz et al., 2006). The other mechanism is hypoxic resting syndrome, where the glycolysis and fermentation processes are accelerated to maintain the energy level, while the tricarboxylic acid cycle (TCA) is inhibited (Bailey-Serres and Voesenek, 2008; Mustroph et al., 2010; Bailey-Serres et al., 2012; Pan et al., 2021). Several life processes, such as cell differentiation and growth, protein synthesis, and cell wall formation, are also slowed down. Plants can survive this energy crisis by reducing energy consumption (Voesenek and Bailey-Serres, 2015). In our previous study, the morphological characteristics of P. ostii did not change under waterlogging stress (Hu et al., 2017). In the present study, the genes encoding ADH, GAPDH, PDC, and PK involved in glycolysis and fermentation were significantly upregulated under waterlogging as identified by RNA-Seq and qRT-PCR verification.
We also found that the genes encoding IDH and MDH involved in TCA cycle were significantly downregulated under waterlogging stress. Moreover, most DEGs involved in carotenoid biosynthesis, monoterpenoid biosynthesis, starch and sucrose metabolism, phenylpropanoid biosynthesis, and isoflavonoid biosynthesis were downregulated. Our results supported the longstanding notion that plants respond to waterlogging stress by regulating energy production and consumption (Bailey-Serres et al., 2012). As for the waterlogging recovery treatment, the altered gene expression trends mentioned above were opposite to waterlogging treatment. All these results collectively indicated that the strategy of P. ostii adopted hypoxic resting syndrome in response to WT.
Hypoxia caused by waterlogging stress modulates the carbon metabolism involved in energy production, consumption, and nitrogen metabolism (Ren et al., 2017). Nitrate reductase and nitrite reductase are two critical enzymes involved in the regulation metabolism, and they can convert nitrate in plants to nitric oxide (Juntawong et al., 2014). However, in the present study, the genes of these two enzymes were not differentially regulated in WT and WRT groups. We also found that the gene encoding nitrate transporter was significantly downregulated under WT and upregulated in response to WRT, as verified by qRT-PCR. Nitrate transporter is responsible for nitrate uptake by plant roots from the soil, its transport, and intracellular redistribution in plants (Jia et al., 2014). We hypothesized that the roots of P. ostii, under waterlogging stress, reduce the production of endogenous nitric oxide by decreasing nitrate uptake from the soil, which is important for avoiding the damage caused by a higher concentration of nitric oxide in the plants.
Aquaporins, proteins belonging to Membrane Intrinsic Proteins (MIP) family facilitate the bidirectional transport of water through biological membranes (Martinez-Ballesta et al., 2014). Enabling enhanced water absorption of roots is one of the main mechanisms for maintaining water content under adverse conditions. Aquaporins play a major role in regulating the hydraulic conductivity that ultimately affects the water uptake capacity of plants (Chaumont and Tyerman, 2014). Numerous studies have been conducted to elucidate the important roles of aquaporins in several plant species in response to salt and drought stresses (Kadam et al., 2017), while there is a severe lack of adequate attention to aquaporins in response to waterlogging stresses, especially for Paeonia species. The gene encoding TIP1-3 was significantly downregulated in the WT group, followed by a significant increase after WRT, although it did not return to the level of the CK group. In sorghum roots, the observed initial upregulation of aquaporin gene expression in response to short-term exposure (18 h) to stress may enhance water uptake to maintain the plant water status, while reduced expression after prolonged exposure (96 h) to stress may reduce hydraulic conductivity (Lv et al., 2016). In this study, P. ostii suffered a long-term waterlogging (3 d), and the water content in the roots was nearly saturated. Hence, the aquaporins gene expression naturally decreased to reduce the water intake, which was in line with the above result. After 7 d of waterlogging recovery, the aquaporin gene was significantly upregulated, but it did not reach the level of CK, which could be attributed to P. ostii still being in the recovery phase.
Ethylene, an important plant hormone, mediates tolerance and adaption of plants to adverse abiotic stresses including waterlogging (Lv et al., 2016). In the ethylene biosynthesis pathway, S-adenosylmethionine is used as a substrate to generate 1-aminocyclopropane-1-carboxylate (ACC) under the catalysis of ACC synthase (ACS), and then ACC gets oxidized by ACC oxidase (ACO) to ethylene. Therefore, ACS and ACO are two key enzymes controlling ethylene biosynthesis in plants (Zhang et al., 2018b). Under the WT condition, the expression of ACO in the root of P. ostii was significantly upregulated as verified by qRT-PCR, and downregulated in response to the WRT condition. Unexpectedly, ACS was not a DEG in response to WT and WRT conditions, which may be because the production of ACC did not take place in the roots of P. ostii. In contrast to our results, ACS expression was significantly upregulated, while ACO was not expressed differentially in the waterlogged Rehmannia glutinosa roots (Wang et al., 2017). This observation indicated that the oxidation of ACC to ethylene was not in the roots of R. glutinosa, and hence the difference of ethylene production in response to waterlogging stress between different plant species.
Transcription factors (TFs) are master regulators of abiotic stress responses in plants. These TFs probably initiate the indirect-late phase of responses by binding to cis-acting elements in the promoters of specific target genes encoding proteins with specific functions (Lindemose et al., 2013). In the present study, 3.80% and 4.40% of the DEGs encoding TFs were identified in the two treatments respectively, indicating that transcriptional regulation played an important role in response to WT and WRT conditions. Consistent with the previous finding that ethylene synthesis and perception were activated (Juntawong et al., 2014), our study recognized AP2/EREBP as the most pronounced TFs in response to both WT and WRT. Their expression patterns were closely followed by the MYB TFs, which were regarded as the second largest TF family in response to WT and WRT. Some studies have reported that MYB TFs played a critical regulatory role in response to waterlogging stress of many plants (Wang et al., 2017; Li et al., 2018) and various stress of P. ostii (Gai et al., 2013; Wang et al., 2016). bHLH and WRKY TFs were expressed more in response to WT, but less in WRT, which suggested that these two TFs may play a more critical role in waterlogging stress. Four TFs genes encoding AP2/EREBP, WRKY, MYB, and NAC were chosen for further expression analysis by qRT-PCR due to their high differential expression as captured in the sequencing results. They were all induced in response to the WT condition, while the decreasing expression levels were observed under the WRT condition. These results indicated that all the validated TFs were highly responsive to waterlogging and waterlogging recovery. These TFs need to be analyzed further to identify their target genes, which may thus provide important insights into the molecular mechanisms underlying waterlogging tolerance in P. ostii.
In summary, this study was the first report on a comprehensive physiological and transcriptomic analysis of P. ostii roots in response to waterlogging and its recovery. The strategy of P. ostii adopted hypoxic resting syndrome in response to waterlogging. A hybrid approach combining SMRT and NGS sequencing platforms was applied in the present study, and the full-length transcriptome sequencing of P. ostii was reported for the first time. A lot of responsive transcripts and genes that might play an important role in waterlogging and its recovery of P. ostii were identified, especially those involved in the glycolysis and fermentation, tricarboxylic acid cycle, nitrogen metabolism, water absorption, ethylene biosynthesis, and transcriptional regulation. The detailed characterization of individual genes must be performed as the next step to clarify their specific functions. These results will enhance our understanding of defense mechanisms against waterlogging stress in P. ostii, and the full-length transcriptome established in the present study will provide the data support for future studies on P. ostii at the molecular levels.
Data availability statement
The datasets presented in this study can be found in online repositories. The name of the repository and accession number are as follows: Genome Sequence Archive, CRA004511.
Author contributions
JY and YH conceived the project. JY designed all the experiments. JY, XZ, and XL analyzed the data. XZ and MZ collected the samples and performed the experiments. XZ, JY, and YH drafted and revised the manuscript. All authors read and approved the final version of the manuscript.
Funding
This research was funded by Shanghai Science and Technology Committee (14JC1403902, 21DZ1202000), Special Fund for Scientific Research of Shanghai Landscaping Administration Bureau Program (G222415), National Natural Science Foundation of China (31470328, 32000240), Open Project of Shanghai Key Laboratory of Plant Functional Genomics and Resources (PFGR202203), and Key project at the central government level: The ability establishment of sustainable use for valuable Chinese medicine resources (2060302).
Acknowledgments
We thank the joint Postgraduate Shaobo Du, Yanli Hu, Yunpeng Wang, and Yanan Xuan for collecting samples.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.1030584/full#supplementary-material
References
Altschul, S. F., Gish, W., Miller, W., Myers, E. W., Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2
Bailey-Serres, J., Fukao, T., Gibbs, D. J., Holdsworth, M. J., Lee, S. C., Licausi, F., et al. (2012). Making sense of low oxygen sensing. Trends Plant Sci. 17, 129–138. doi: 10.1016/j.tplants.2011.12.004
Bailey-Serres, J., Voesenek, L. A. C. J. (2008). Flooding stress: acclimations and genetic diversity. Annu. Rev. Plant Biol. 59, 313–339. doi: 10.1146/annurev.arplant.59.032607.092752
Chang, Y., Hu, T., Zhang, W., Zhou, L., Wang, Y., Jiang, Z. (2019). Transcriptome profiling for floral development in reblooming cultivar ‘High noon’ of Paeonia suffruticosa. Sci. Data 6, 217. doi: 10.1038/s41597-019-0240-1
Chao, Q., Gao, Z.-F., Zhang, D., Zhao, B.-G., Dong, F.-Q., Fu, C.-X., et al. (2019). The developmental dynamics of the Populus stem transcriptome. Plant Biotechnol. J. 17, 206–219. doi: 10.1111/pbi.12958
Chaumont, F., Tyerman, S. D. (2014). Aquaporins: highly regulated channels controlling plant water relations. Plant Physiol. 164, 1600–1618. doi: 10.1104/pp.113.233791
Chen, W., Yao, Q., Patil, G. B., Agarwal, G., Deshmukh, R. K., Lin, L., et al. (2016). Identification and comparative analysis of differential gene expression in soybean leaf tissue under drought and flooding stress revealed by RNA-seq. Front. Plant Sci. 7 1044. doi: 10.3389/fpls.2016.01044
Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Deng, Y., Zheng, H., Yan, Z., Liao, D., Li, C., Zhou, J., et al. (2018). Full-length transcriptome survey and expression analysis of Cassia obtusifolia to discover putative genes related to aurantio-obtusin biosynthesis, seed formation and development, and stress response. Int. J. Mol. Sci. 19, E2476. doi: 10.3390/ijms19092476
Ding, J., Liang, P., Wu, P., Zhu, M., Li, C., Zhu, X., et al. (2020). Effects of waterlogging on grain yield and associated traits of historic wheat cultivars in the middle and lower reaches of the Yangtze river, China. Field Crop Res. 246, 107695. doi: 10.1016/j.fcr.2019.107695
Fu, L., Niu, B., Zhu, Z., Wu, S., Li, W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565
Gai, S., Zhang, Y., Liu, C., Zhang, Y., Zheng, G. (2013). Transcript profiling of Paoenia ostii during artificial chilling induced dormancy release identifies activation of GA pathway and carbohydrate metabolism. PloS One 8, e55297. doi: 10.1371/journal.pone.0055297
Glenz, C., Schlaepfer, R., Iorgulescu, I., Kienast, F. (2006). Flooding tolerance of central European tree and shrub species. For. Ecol. Manag 235, 1–13. doi: 10.1016/j.foreco.2006.05.065
Hu, Y., Han, J. (2018). Tree peony, its varieties, cultivation and use in East China (Beijing: Science Press).
Hu, Y., Yuan, J., Hu, Y. (2017). Effects of waterlogging and recovery on the photosynthetic characteristics and root activity of Paeonia ostii. Acta Bot. Boreal-Occident Sin. 37, 1586–1594. doi: 10.7606/j.issn.10004025.2017.08.1586
Hu, Z., Zhang, Y., He, Y., Cao, Q., Zhang, T., Lou, L., et al. (2020). Full-length transcriptome assembly of Italian ryegrass root integrated with RNA-seq to identify genes in response to plant cadmium stress. Int. J. Mol. Sci. 21, E1067. doi: 10.3390/ijms21031067
Jia, D., Wang, Y., Liu, Y., Hu, J., Guo, Y., Gao, L., et al. (2018). SMRT sequencing of full-length transcriptome of flea beetle Agasicles hygrophila (Selman and vogt). Sci. Rep. 8, 2197. doi: 10.1038/s41598-018-20181-y
Jia, H., Zhang, H., Liu, W., Cui, H., Liu, G. (2014). Function and regulation mechanisms of nitrate transporters in higher plants. Biotech. Bull. 6, 14–21. doi: 10.13560/j.cnki.biotech.bull.1985.2014.06.003
Jones, P., Binns, D., Chang, H.-Y., Fraser, M., Li, W., McAnulla, C., et al. (2014). InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240. doi: 10.1093/bioinformatics/btu031
Juntawong, P., Sirikhachornkit, A., Pimjan, R., Sonthirod, C., Sangsrakru, D., Yoocha, T., et al. (2014). Elucidation of the molecular responses to waterlogging in jatropha roots by transcriptome profiling. Front. Plant Sci. 5. doi: 10.3389/fpls.2014.00658
Kadam, S., Abril, A., Dhanapal, A. P., Koester, R. P., Vermerris, W., Jose, S., et al. (2017). Characterization and regulation of aquaporin genes of sorghum [Sorghum bicolor (L.) moench] in response to waterlogging stress. Front. Plant Sci. 8. doi: 10.3389/fpls.2017.00862
Kim, D., Langmead, B., Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317
Langan, P., Bernád, V., Walsh, J., Henchy, J., Khodaeiaminjan, M., Mangina, E., et al. (2022). Phenotyping for waterlogging tolerance in crops: current trends and future prospects. J. Exp. Bot. 73 (15), 5149–5169. doi: 10.1093/jxb/erac243
Langmead, B., Salzberg, S. L. (2012). Fast gapped-read alignment with bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Lee, Y.-H., Kim, K.-S., Jang, Y.-S., Hwang, J.-H., Lee, D.-H., Choi, I.-H. (2014). Global gene expression responses to waterlogging in leaves of rape seedlings. Plant Cell Rep. 33, 289–299. doi: 10.1007/s00299-013-1529-8
Lefevre, M., Kris-Etherton, P. M., Zhao, G., Tracy, R. P. (2004). Dietary fatty acids, hemostasis, and cardiovascular disease risk. J. Am. Diet Assoc. 104, 410–419. doi: 10.1016/j.jada.2003.12.022
Li, B., Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinf. 12, 323. doi: 10.1186/1471-2105-12-323
Li, R., Li, Y., Kristiansen, K., Wang, J. (2008). SOAP: short oligonucleotide alignment program. Bioinformatics 24, 713–714. doi: 10.1093/bioinformatics/btn025
Lindemose, S., O’Shea, C., Jensen, M. K., Skriver, K. (2013). Structure, function and networks of transcription factors involved in abiotic stress responses. Int. J. Mol. Sci. 14, 5842–5878. doi: 10.3390/ijms14035842
Li, Y., Yang, J., Liu, R., Zhou, D., Gan, L., Wu, Y. (2018). Transcriptome analysis of transcription factor of suitable Cynodon dactylon in three gorges reservoir area under waterlogging. Southwest China J. Agr Sci. 31, 265–269. doi: 10.16213/j.cnki.scjas.2018.2.008
Li, S., Yuan, R., Chen, L., Wang, L., Hao, X., Wang, L.-J., et al. (2015). Systematic qualitative and quantitative assessment of fatty acids in the seeds of 60 tree peony (Paeonia section Moutan DC.) cultivars by GC–MS. Food Chem. 173, 133–140. doi: 10.1016/j.foodchem.2014.10.017
Li, J., Zhang, X., Zhao, X. (2011). Tree peony of China (Beijing: Encyclopedia of China Publishing House).
Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi: 10.1186/s13059-014-0550-8
Lv, Y., Fu, S., Chen, S., Zhang, W., Qi, C. (2016). Ethylene response factor BnERF2-like (ERF2.4) from brassica napus l. enhances submergence tolerance and alleviates oxidative damage caused by submergence in Arabidopsis thaliana. Crop J. 4, 199–211. doi: 10.1016/j.cj.2016.01.004
Martinez-Ballesta, M., del, C., Carvajal, M. (2014). New challenges in plant aquaporin biotechnology. Plant Sci. 217–218, 71–77. doi: 10.1016/j.plantsci.2013.12.006
McDonald, M. P., Galwey, N. W., Colmer, T. D. (2002). Similarity and diversity in adventitious root anatomy as related to root aeration among a range of wetland and dryland grass species. Plant Cell Environ. 25, 441–451. doi: 10.1046/j.0016-8025.2001.00817.x
Minoche, A. E., Dohm, J. C., Schneider, J., Holtgräwe, D., Viehöver, P., Montfort, M., et al. (2015). Exploiting single-molecule transcript sequencing for eukaryotic gene prediction. Genome Biol. 16, 184. doi: 10.1186/s13059-015-0729-7
Mustroph, A., Lee, S. C., Oosumi, T., Zanetti, M. E., Yang, H., Ma, K., et al. (2010). Cross-kingdom comparison of transcriptomic adjustments to low-oxygen stress highlights conserved and plant-specific responses. Plant Physiol. 152 (3), 1484–1500. doi: 10.1104/pp.109.151845
Pang, T., Ye, C.-Y., Xia, X., Yin, W. (2013). De novo sequencing and transcriptome analysis of the desert shrub, ammopiptanthus mongolicus, during cold acclimation using Illumina/Solexa. BMC Genomics 14, 488. doi: 10.1186/1471-2164-14-488
Pan, J., Sharif, R., Xu, X., Chen, X. (2021). Mechanisms of waterlogging tolerance in plants: Research progress and prospects. Front. Plant Sci. 11. doi: 10.3389/fpls.2020.627331
Pedersen, O., Sauter, M., Colmer, T. D., Nakazono, M. (2021). Regulation of root adaptive anatomical and morphological traits during low soil oxygen. New Phytol. 229, 42–49. doi: 10.1111/nph.16375
Postnikova, O. A., Shao, J., Nemchinov, L. G. (2013). Analysis of the alfalfa root transcriptome in response to salinity stress. Plant Cell Physiol. 54, 1041–1055. doi: 10.1093/pcp/pct056
Ren, C.-G., Kong, C.-C., Yan, K., Zhang, H., Luo, Y.-M., Xie, Z.-H. (2017). Elucidation of the molecular responses to waterlogging in Sesbania cannabina roots by transcriptome profiling. Sci. Rep. 7, 9256. doi: 10.1038/s41598-017-07740-5
Shahidi, F., Miraliakbari, H. (2005). Omega-3 fatty acids in health and disease: Part 2–health effects of omega-3 fatty acids in autoimmune diseases, mental health, and gene expression. J. Med. Food 8, 133–148. doi: 10.1089/jmf.2005.8.133
Sharon, D., Tilgner, H., Grubert, F., Snyder, M. (2013). A single-molecule long-read survey of the human transcriptome. Nat. Biotechnol. 31, 1009–1014. doi: 10.1038/nbt.2705
Steijger, T., Abril, J. F., Engström, P. G., Kokocinski, F., Hubbard, T. J., Guigó, R., et al. (2013). Assessment of transcript reconstruction methods for RNA-seq. Nat. Methods 10, 1177–1184. doi: 10.1038/nmeth.2714
Voesenek, L. A. C. J., Bailey-Serres, J. (2015). Flood adaptive traits and processes: an overview. New Phytol. 206, 57–73. doi: 10.1111/nph.13209
Wang, Y., Dong, C., Xue, Z., Jin, Q., Xu, Y. (2016). De novo transcriptome sequencing and discovery of genes related to copper tolerance in Paeonia ostii. Gene 576, 126–135. doi: 10.1016/j.gene.2015.09.077
Wang, C., Li, X., Wang, X., Li, M., Zhang, Z., Chen, X. (2017). Identification of key genes in response to waterlogging stress in root of Rehmannia glutinosa. Mod Chin. Med. 19, 232–238. doi: 10.13313/j.issn.1673-4890.2017.2.015
Wu, S., Jin, X., Zhang, M., Zhang, F., Luo, F. (2009). Effects of exogenous salicylic acid on heat tolerance of tree peony seedlings under high temperature stress. J. Henan Agr. Sci. 47, 98–103. doi: 10.15933/j.cnki.1004-3268.2018.06.017
Yeung, E., van Veen, H., Vashisht, D., Sobral Paiva, A. L., Hummel, M., Rankenberg, T., et al. (2018). A stress recovery signaling network for enhanced flooding tolerance in Arabidopsis thaliana. Proc. Natl. Acad. Sci. U.S.A. 115, E6085–E6094. doi: 10.1073/pnas.1803841115
Yin, D., Sun, D., Han, Z., Ni, D., Norris, A., Jiang, C.-Z. (2019). PhERF2, an ethylene-responsive element binding factor, plays an essential role in waterlogging tolerance of petunia. Hortic. Res. 6, 83. doi: 10.1038/s41438-019-0165-z
Yu, S., Du, S., Yuan, J., Hu, Y. (2016). Fatty acid profile in the seeds and seed tissues of Paeonia l. species as new oil plant resources. Sci. Rep. 6, 26944. doi: 10.1038/srep26944
Zhang, Z., Luan, D., Chen, G., Xiao, F., Yin, D., Ni, D., et al. (2018b). Isolation and characterization of a waterlogging-responsive gene involved in ethylene biosynthesis in chrysanthemum. Acta Physiol. Plant 40, 100. doi: 10.1007/s11738-018-2667-y
Zhang, X., Shi, Q., Ji, D., Niu, L., Zhang, Y. (2017a). Determination of the phenolic content, profile, and antioxidant activity of seeds from nine tree peony (Paeonia section Moutan DC.) species native to China. Food Res. Int. 97, 141–148. doi: 10.1016/j.foodres.2017.03.018
Zhang, X., Zhang, G., Jin, M., Niu, L., Zhang, Y. (2018a). Variation in phenolic content, profile, and antioxidant activity of seeds among different Paeonia ostii cultivated populations in China. Chem. Biodivers 15, e1800093. doi: 10.1002/cbdv.201800093
Zhang, X., Zhang, Y., Niu, L., Sun, J., Li, L., Zhang, J., et al. (2017b). Chemometric classification of different tree peony species native to china based on the assessment of major fatty acids of seed oil and phenotypic characteristics of the seeds. Chem. Biodivers 14, e1600111. doi: 10.1002/cbdv.201600111
Keywords: tree peony, transcriptome, waterlogging, waterlogging recovery, hypoxia
Citation: Zhang X, Liu X, Zhou M, Hu Y and Yuan J (2022) PacBio full-length sequencing integrated with RNA-seq reveals the molecular mechanism of waterlogging and its recovery in Paeonia ostii. Front. Plant Sci. 13:1030584. doi: 10.3389/fpls.2022.1030584
Received: 29 August 2022; Accepted: 12 October 2022;
Published: 03 November 2022.
Edited by:
Tangchun Zheng, Beijing Forestry University, ChinaReviewed by:
Minhuan Zhang, Central South University Forestry and Technology, ChinaShunli Wang, Chinese Academy of Agricultural Sciences (CAAS), China
Donghui Peng, Fujian Agriculture and Forestry University, China
Copyright © 2022 Zhang, Liu, Zhou, Hu and Yuan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Junhui Yuan, eXVhbmp1bmh1aWdzbHlAMTI2LmNvbQ==; Yonghong Hu, aHV5b25naG9uZ0Bjc25iZ3NoLmNu
†These authors have contributed equally to this work