De novo sequencing, assembly, and analysis of the Taxodium‘Zhongshansa’ roots and shoots transcriptome in response to short-term waterlogging

Abstract

Background

Taxodium is renowned for its strong tolerance to waterlogging stress, thus it has great ecological and economic potential. However, the scant genomic resources in genus Taxodium have greatly hindered further exploration of its underlying flood-tolerance mechanism. Taxodium ‘Zhongshansa’ is an interspecies hybrid of T. distichum and T. mucronatum, and has been widely planted in southeastern China. To understand the genetic basis of its flood tolerance, we analyzed the transcriptomes of Taxodium ‘Zhongshansa’ roots and shoots in response to short-term waterlogging.

Results

RNA-seq was used to analyze genome-wide transcriptome changes of Taxodium ‘Zhongshansa 406’ clone root and shoot treated with 1 h of soil-waterlogging stress. After de novo assembly, 108,692 unigenes were achieved, and 70,260 (64.64%) of them were annotated. There were 2090 differentially expressed genes (DEGs) found in roots and 394 in shoots, with 174 shared by both of them, indicating that the aerial parts were also affected. Under waterlogging stress, the primary reaction of hypoxic-treated root was to activate the antioxidative defense system to prevent cells experiencing reactive oxygen species (ROS) poisoning. As respiration was inhibited and ATP decreased, another quick coping mechanism was repressing the energy-consuming biosynthetic processes through the whole plant. The glycolysis and fermentation pathway was activated to maintain ATP production in the hypoxic root. Constantly, the demand for carbohydrates increased, and carbohydrate metabolism were accumulated in the root as well as the shoot, possibly indicating that systemic communications between waterlogged and non-waterlogged tissues facilated survival. Amino acid metabolism was also greatly influenced, with down-regulation of genes involvedin serine degradation and up-regulation of aspartic acid degradation. Additionally, a non-symbiotic hemoglobin class 1 gene was up-regulated, which may also help the ATP production. Moreover, the gene expression pattern of 5 unigenes involving in the glycolysis pathway revealed by qRT-PCR confirmed the RNA-Seq data.

Conclusions

We conclude that ROS detoxification and energy maintenance were the primary coping mechanisms of ‘Zhongshansa’ in surviving oxygen deficiency, which may be responsible for its remarkable waterlogging tolerance. Our study not only provided the first large-scale assessment of genomic resources of Taxodium but also guidelines for probing the molecular mechanism underlying ‘Zhongshansa’ waterlogging tolerance.

Background

The genus Taxodium is historically recognized as containing three species: T. distichum (baldcypress), T. mucronatum (Montezuma cypress) and T. ascendens (pondcypress) [1]. However, there is still some debate concerning the taxonomy of these three taxa [1]. In the present study, we take the taxonomic opinion of Zheng [2] who treated the genus Taxodium as three distinct species. Taxodium are extremely flood-tolerant conifers in the cypress family, and thus have many positive environmental attributes both as wetland species [3] and as landscape plants [1].

To develop optimal woody plants for afforestation in the coastal and wetland areas of southeastern China, a number of interspecies crosses among the three Taxodium species have been conducted since the 1970s, from which a batch of superior hybrid clones have been selected, such as ‘Zhongshansa 302’ (T. distichum × T. mucronatum), ‘Zhongshansa 118’ [(T. distichum × T. mucronatum) × T. mucronatum] and ‘Zhongshansa 406’ (T. mucronatum × T. distichum) [4]. Taxodium ‘Zhongshansa’ are conical, deciduous to semi-evergreen conifers with needle-like leaves, and are interspecies hybrids of T. mucronatum and T. distichum. ‘Zhongshansa’ are extremely tolerant to waterlogging [4] and can survive for months with their roots in flooded soil where most tree species cannot subsist. Currently in southeastern China, ‘Zhongshansa’ have been widely used as timber trees in river network areas, as windbreak trees in coastal areas and as landscape trees in urban areas. Despite its great ecological and economic potential, genomic information on genus Taxodium is scarce, which greatly hinders the development of molecular markers, further exploration of its underlying flood-tolerance mechanism and other genetic research.

Higher plants are aerobic organisms. Since the diffusion rate of molecular oxygen in water is much lower than in air, soil waterlogging is a serious obstacle to plant growth and development, which may make plants hypoxic or anoxic. The response of plants to external hypoxia has been intensively studied in the past. Proteomics research has identified a set of about 20 anaerobically induced polypeptides (ANPs) [5]. ANPs have been demonstrated as essential for tolerance to low oxygen in a number of plant species [6],[7]. Further studies showed that the majority of ANPs were involved in the glycolysis and fermentation pathways [8]. Subsequently, microarray studies have been performed on the low-oxygen response in Arabidopsis thaliana[9], maize [10], cotton [5], poplar [11] and other plants. All of these rapid changes in a large number of transcripts involving not only well-known ANPs [12], but also those previously unknown to be involved in hypoxia or anoxia response, indicating that plants have complex responses to low oxygen [5],[13].

Compared with microarrays, the RNA-Seq approach has higher sensitivity which includes both low- and high-level gene expression [14]. These advantages have resulted in the increased application of RNA-seq to elucidate the response of plants to various environmental stresses, such as cold [15], salt [16],[17] and drought [16],[18]. RNA-seq has also been successfully used in crops’ responses to waterlogging stress, such as maize [19], cucumber [20], sesame [21] and rape [22]. However, seldom reports has been found on the woody plants.

To better understand the molecular mechanisms of the response of ‘Zhongshansa’ to soil waterlogging, the global gene transcription changes in both submerged roots and aerial shoot tissues of waterlogged Taxodium ‘Zhongshansa 406’ clone were examined using the Illumina HiSeq™ 2000 sequencing platform (Illumina Inc., San Diego, CA, USA). We focused on the early stage of ‘Zhongshansa’ response to waterlogging stress because it determines the switch from normal to low-oxygen metabolism and plays an essential role in plant survival [8]. To our knowledge, this is the first large-scale assessment of Taxodium genomic resources. Our results will facilitate understanding of the response of flood-tolerant woody plants to soil waterlogging stress.

Methods

Plant growth and water treatments

Cuttings of the Taxodium clone ‘Zhongshansa 406’ were cultured in plastic pots in a ventilated greenhouse of the Nanjing Botanical Garden in April 2010. In July 2013, six plantlets were moved and cultured at room temperature (approximately 20°C), using a photoperiod of 16/8 h of light/dark. Two weeks later, plantlets were divided evenly into two groups: one served as the control sample (CK), while the other was treated with tap water with the plastic pots immersed as the waterlogging treated sample (CT). The roots and shoots of CT were sampled at 1 h after the application of fresh water. The roots and shoots of CK were also sampled at the same time-point. The primary root with some lateral roots and the shoot apex with three leaves were simultaneously collected from each individual plant, and were separately frozen in liquid nitrogen and stored at −80°C prior to RNA extraction. Roots were washed carefully to prevent mechanical damage. In total, 4 RNA pools were achieved, e.g. CK root, CK shoot, CT root and CT shoot, and each of the RNA pool was made by the mixture of the same tissues from 3 plantlets in the same group.

RNA isolation, cDNA library construction and sequencing

Total RNA of roots was first crudely extracted using the RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (Tiangen, Beijing, China), and then purified with the RNA Clean-up Kit (Tiangen). For the leaves, total RNA was isolated with the PLANTeasy Plant RNA Extraction Kit (Yuanpinghao, Beijing, China) according to the manufacturer’s instructions. RNA quality detection, cDNA library construction and Illumina deep sequencing were performed following the method of Lv [23], and 150 bp paired-end reads were generated.

Assembly and annotation

To get high-quality clean reads, in-house perl scripts were used to process raw data, which removed reads containing adapters, low-quality reads and reads containing poly-N. The calculation of Q20, Q30, GC-content and sequence duplication level, and other downstream analyses were based on the clean reads. Transcriptome assembly was achieved using Trinity [24].

Gene function was annotated based on the following seven databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database) and GO (Gene Ontology), using BLAST with a cutoff E-value of 10⁻⁵.

Quantification of gene expression levels and differential expression analysis

Gene expression levels were estimated by RSEM [25] for each sample. Clean data were mapped back onto the assembled transcriptome. Readcount for each gene was obtained from the mapping results and normalized to reads per kb of exon model per million mapped reads (RPKM). Prior to differential gene expression analysis for each sequenced library, the readcounts were adjusted by edgeR program package [26] through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq (2010) R package. P-value was adjusted using q-value [27]; with q-value < 0.005 and |log2 (foldchange)| > 1 as the threshold for significant differential expression. GO enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GO seq R packages based Wallenius non-central hyper-geometric distribution [28], which can adjust for gene length bias in DEGs. KOBAS [29] software was used to test the statistical enrichment of DEGs in KEGG pathways.

qRT-PCR analysis

The expression patterns of five genes involving in the glycolysis pathway (Gene ID: comp71558_c0, comp63755_c0, comp75584_c0, comp53892_c1, and comp62913_c0) were analyzed using qRT-PCR. New plant materials of the same clone were used for the RNA extraction for the qRT-PCR assays. And three biological replicates were made. Gene-specific primers were designed according to the reference unigene sequences using the Primer Premier 5.0. A HiScript^TM Q RT SuperMix for qPCR (Vazyme, Nanjing, China) was used to synthesize the cDNAs and real-time quantification was performed using a ABI StepOne^TM Plus system and the AceQ^TM qPCR SYBR® Green Master Mix (Vazyme, Nanjing, China). PCR cycling was denatured using a program of 95°C for 5 min, and 40 cycles of 95°C for 10 s and 60°C for 30 s ‘Zhongshansa’ actin gene (forward: 5′- TTAACATTGTGACCTGTGCGAACT - 3′, and reverse: 5′-ACAACAAGGAAAGTATAGCCAGCAA −3′) was used as a normalizer, and the relative expression levels of genes were presented by 2^-△△CT as all the genes tested show highly similar amplification efficiency around 0.95 (Additional file 1).

Results

Transcriptome sequencing and assembly

Illumina sequencing data from ‘Zhongshansa’ roots and shoots were deposited in the NCBI SRA database under accession number SRP043177. In total, 174,958,744 Illumina PE raw reads were generated (Table 1). After removing adaptor sequences, ambiguous nucleotides and low-quality sequences, there were 153,993,822 million clean reads remaining. Assembly of clean reads resulted in 108,692 unigenes in the range of 201–14,489 bp with a N50 length of 1123 bp (Figure 1).

Table 1 Summary of sequences analysis

Figure 1

Length distribution of assembled unigenes.

Sequence annotation

The unigenes were annotated by aligning with the seven public databases (Table 2). Analyses showed that 61,087 unigenes (56.2%) had significant matches in the Nr database, 21,203 (19.5%) in the Nt database and 44,761 (41.18%) in the Swiss-Prot database. In total, there were 70,260 unigenes (64.64%) successfully annotated in at least one of the Nr, Nt, Swiss-Prot, KEGG, GO, COG and Pfam databases, with 7622 unigenes (7.01%) in all seven databases.

Table 2 BLAST analysis of non-redundant unigenes against public databases

For GO analysis, there were 50,929 unigenes divided into three ontologies (Figure 2). For biological process (BP) category, genes involved in ‘cellular process’ (28,970), ‘metabolic process’ (28,659) and ‘single-organism process’ (13,853) were highly represented. The cellular component (CC) category mainly comprised proteins involved in ‘cell’ (17,488), ‘cell part’ (17,471) and ‘organelle’ (11,813). Within the molecular function (MF) category, ‘binding’ (28,115), ‘catalytic activity’ (25,271) and ‘transporter activity’ (4135) were highly represented.

Figure 2

GO categorization of non-redundant unigenes.

In addition, all unigenes were subjected to a search against the COG database for functional prediction and classification. In total, there were 31,506 unigenes assigned to COG classification and divided into 26 specific categories (Figure 3). The ‘general functional prediction only’ (4714) was the largest group, followed by ‘post-translational modification, protein turnover, chaperon’ (4440), ‘translation’ (2959), ‘signal transduction’ (2755) and ‘energy production and conversion’ (2074). Only a few unigenes were assigned to ‘extracellular structures’ (177) and ‘cell motility’ (31).

Figure 3

COG annotation of putative proteins.

The unigene metabolic pathway analysis was also conducted using the KEGG annotation system. This process predicted a total of 258 pathways, representing a total of 22,871 unigenes (Figure 4). The pathways involving the highest number of unique transcripts were ‘translation’ (2743), followed by ‘carbohydrate metabolism’ (2646) and ‘energy metabolism’ (2176).

Figure 4

KEGG annotation of putative proteins.

Differential expression analysis of assembled ‘Zhongshansa’ transcripts under waterlogging treatments in different tissues

Differential expression analysis was firstly performed between the two tissues. DEGs (q-value < 0.005 and |log2 (foldchange)| >1) were defined as genes that were significantly enriched or depleted in one tissue relative to the other tissue. In the CK, there were 4730 DEGs between the shoots and roots, and 4677 DEGs between treated shoots and roots.

Then, the DEGs between the CK and CT were analyzed. Of 108,692 (2.1%) unigenes, 2310 were identified as DEGs in at least one tissue between CT and CK plants (Figure 5). Among them, 2090 DEGs were found in roots and 394 in shoots. In this study, DEGs with higher expression levels in CT compared with CK were denoted as ‘up-regulated’, while those with lower expression levels in CT were ‘down-regulated’. There were 174 DEGs shared by both tissues, among which 28 showed opposite trends in expression between roots and shoots, with 10 up-regulated and 18 down-regulated in roots. The remaining 146 DEGs showed similar expression differences in each tissue, including 99 down-regulated and 47 up-regulated DEGs.

Figure 5

Venn diagrams of the differential expression transcripts under waterlogging treatment in root and leaf samples. The numbers of DEGs exclusively up- or down-regulated in one tissue are shown in each circle. The numbers of DEGs with a common or opposite tendency of expression changes between the two tissues are shown in the overlapping regions. The total numbers of up- or down-regulated genes in each tissue are shown outside the circles.

There were 1916 genes exclusively differentially expressed in roots, with 1009 down-regulated and 907 up-regulated. There were 220 DEGs (167 up-regulated and 53 down-regulated) exclusively changed in shoots.

Functional classification of DEGs

To further characterize the expression changes discussed above, we conducted GO enrichment analysis for DEGs with the whole transcriptome as the background. GO analysis was conducted on the DEGs between the shoot and root in CK. GO enrichment analysis of the up-regulated DEGs in the shoot compared to root indicated some shoot-specific or strongly performed functions. mRNAs in the shoot were highly enriched encoding proteins involved in all aspects of photosynthesis, with ‘photosynthesis’, ‘oxidation-reduction process’, ‘photosynthesis, light reaction’ and ‘photosynthesis, light harvesting’ listed as the top-four enriched BPs. Research in Arabidopsis indicated that genes associated with photosynthesis were abundantly expressed in the photosynthetic cells and guard cells of shoots, while largely absent from root mRNAs [30] – this is exactly what our data also suggests. The following highly enriched BPs (corrected p-value <0.005) included processes involved in the biosynthetic and metabolic processesoflipids, steroids, isoprenoids, monocarboxylic acids, fatty acids and others. This indicated that processes involving a series of lipids were quite active in the shoot. Our findings will facilitate research on ‘Zhongshansa’ leaf lipid content.

The down-regulated DEGs in control shoots compared with roots were those whose mRNA were specific or abundant in the root. Not unexpectedly, function categories related to cell proliferation and development were highly enriched, such as ‘negative regulation of growth’, ‘regulation of growth’, ‘cellular component organization or biogenesis’, ‘ribosome biogenesis’, ‘ribonucleoprotein complex biogenesis’ and ‘plant-type cell wall organization’ as the root samples were mainly primary root with some lateral roots. The terms ‘response to stress’, ‘peroxidase reaction’ and ‘response to oxidative stress’ were also among the highly enriched terms. mRNAs were also enriched for binding (heme binding, tetrapyrrole binding, iron ion binding, cation binding and metal ion binding), which is common in the Arabidopsis root [30].

GO analysis was conducted for the up-regulated DEGs in roots (Additional file 2). In the MFcategory, the top three enriched terms were peroxidase activity, oxidoreductase activity acting on peroxide as acceptor, and heme binding. In the CC category, ‘cell wall’, ‘external encapsulating structure’ and ‘apoplast’ were the three dominant enriched terms. In BP, ‘peroxidase reaction’, ‘response to oxidative stress’ and ‘carbohydrate metabolic process’ were the mostly highly enriched. The aspartic metabolism were influenced, with ‘aspartic-type endopeptidase activity’ and ‘aspartic-type peptidase activity’ also enriched (P-value < 0.05). For the down-regulated DEGs in CK compared to CT roots (Additional file 3), ‘ribosome biogenesis’ and ‘ribonucleoprotein complex biogenesis’ were the top-two BPs enriched by the down-regulated DEGs. The ribosome is the place where mRNA istranslated into protein. The decelebration of ribosome and ribonucleoprotein complex biogenesis might imply a great inhibition of protein production in the root. Consistent with this ‘translation’ was the top-four BP enriched by down-regulated DEGs. In treated ‘Zhongshansa’ root, ‘cellular component biogenesis’ and ‘cellular component organization or biogenesis’ were the third and sixth most enriched BPs, respectively. As discussed above, many of the significantly inhibited function categories were highly enriched in control roots. Taken together, proliferation of root cells was greatly limited under hypoxia stress, which may save much energy. The majority of genes involved in mitochondrial electron transport were down-regulated, such as ‘mitochondrial electron transport, cytochrome c to oxygen’, including eight DEGs with seven down-regulated, and the ‘mitochondrial electron transport, NADH to ubiquinone’, with four down-regulated among the six DEGs. Other enriched terms included ‘serine type endopeptidase activity’.

When comparing CT with CK, ‘plant-type cell wall organization’ and ‘plant-type cell wall organization or biogenesis’ were the top-two GO enrichment terms of the down-regulated DEGs in the shoots (Additional file 4) -both of them had four DEGs, which were all repressed. Additionally, 100% of DEGs involvedin ‘cellulose synthase activity’, ‘cellulose synthase (UDP-forming) activity’ and ‘cellulose biosynthetic process’ were also down-regulated. Changes in transcript levels suggested that the energy-demanding cellulose and cell wall biosynthesis processes were greatly inhibited in the shoot. For the GO enrichment analysis of the up-regulated DEGs, ‘transcription, DNA-dependent’, ‘RNA biosynthetic process’ and ‘regulation of gene expression’ were dominant (Additional file 5).

KEGG pathway enrichment analysis for DEGs also revealed both common and tissue-specific patterns of over representations. The top-five enriched pathway by DEGs in CT roots (Additional file 6) (q ≤ 0.05), were phenylpropanoid biosynthesis, phenylalanine metabolism, plant hormone signal transduction, ribosome, protein digestion and absorption. DEGs in shoots were also analyzed (Additional file 7) (q ≤ 0.05). Starch and sucrose metabolism were the top-four enriched pathways by DEGs in CT shoots, compared with CK. There were 468 genes annotated as involved in this pathway, with seven having changed expression under the stress in the shoot. There were two DEGs annotated as encoding trehalose 6-phosphate synthase (TPS)-comp62470_c0 and comp68953_c0-with 5.81- and 2.13-fold increased expression, respectively. Comp64972_c0 encoding a sucrose synthase was 2.29-fold down-regulated in the shoot, which may lead to a slowing of starch production in the root. Four pathways were enriched by DEGs in both tissues: plant hormone signal transduction, carotenoid biosynthesis, starch and sucrose metabolism, and phenylpropanoid biosynthesis.

Perturbation in glycolysisis considered to be the basic characteristic of plant adaption to an aerobic stress [31]. There were 591 unigenes annotated as encoding enzymes involved in glycolysis/gluconeogenesis pathway), with 14 of them differentially expressed between treated and control roots (Figure 6). Most of the DEGs were up-regulated in treated roots,. Two DEGs were annotated as encoding glyceraldehyde 3-phosphate dehydrogenase (GAPDH)-comp68689_c0 and comp64678_c0. And comp64678_c0 was the only down-regulated DEG, indicating that an additional GAPDH isoform may be inhibited in hypoxic root. The activity of the responsible enzyme lactate dehydrogenase (LDH) was up-regulated. Consistent with our results, lactic acid fermentation is activated in the initial stages of root hypoxia in many plants. However, in contrast to animals, the anaerobic metabolism of pyruvate in plants is not limited to the formation of lactate. In gray poplar, LDH transcripts were also rather abundant as an initial reaction to O₂ deprivation, but dropped after about 5 h due to the decrease in cytosolic pH caused by lactic acid [11]. Rather, ethanol is the major fermentation end product for plants. So, the lactic acid fermentation in plant is followed by alcoholic fermentation, with two critical enzymes involved in this process: pyruvate decarboxylase (PDC), which converts pyruvate to acetaldehyde; and alcohol dehydrogenase (ADH), which further metabolizes acetaldehyde to ethanol. In our results, both PDC (comp75584_c0) and ADH (comp71294_c0) were up-regulated in the CT root. By activating alcoholic fermentation, energy was produced in waterlogged ‘Zhongshansa’ root. None of these DEGs showed changed expression in treated shoots. This was consistent with findings for gray poplar [11].

Figure 6

Unigenes predicted to be involved in the glycolysis pathway. Red indicates significantly increased expression in CT compared with CK; green indicates significantly decreased expression; yellow indicates proteins encoded by both up-and down-regulated genes.

Verification of RNA-Seq data by real-time quantitative RT-PCR

To confirm the reliability of the RNA-Seq data, the transcriptional level of 5 unigenes were examined by real-time quantitative PCR (Figure 7). Since, new plant materials were used for the RNA extraction, the fold change did not exactly match the number revealed by the DEG analysis for these genes. All the 5 genes exhibited > 2 fold higher expression in the root in response to waterlogging, while none of them have > 2 fold changes in the shoot. comp53892_c1 annotated as encoding aldehyde dehydrogenase can not be detected in the shoot due to no/low expression, so as the result by the Illumina sequencing technology. Taken together, all the unigenes showed consistent expression patterns that were consistent with the RNA-Seq data, indicating that our experimental results were valid.

Figure 7

Real-time PCR validations of 5 genes in ‘Zhongshsansa’ roots and shoots. Comp53892_c1 was annotated as ALDH, comp63755_c0 was annotated as PFK3, comp71558_c0 was annotated as LDH, comp75584_c0 was annotated as PDC, comp62913_c0 was annotated as PFK2.

Discussion

In this paper, transcriptomes of ‘Zhongshansa 406’ clone roots and shoots were sequenced using the Illumina platform. In total, about 154 million high-quality reads with 23.1 Gb sequence coverage were obtained; there were 108,692 unigenes (≥200 bp) assembled and 64.64% were annotated. As far as we know, this is the first large-scale assessment of Taxodium genomic resources. Our results lay the foundation for development of molecular markers, construction of a genetic map and much other genomics research in Taxodium.

Comparisons of transcriptomes between roots and shoots

We compared the transcriptome differences between root and shoot in the CK. As expected, compared with the root, photosynthesis-relevant mRNAs were abundant in the shoot. The biosynthetic and metabolic processes of a series of lipids were also among the highly enriched, because the leaves of gymnosperms always contain high levels of lipids. Profiling translatomes of discrete cell populations in Arabidopsis showed that all five clusters (clusters 3, 19, 25, 45 and 55) containing the terms ‘response to stress’, ‘peroxidase reaction’ or ‘response to oxidative stress’ – especially cluster 45 was enriched in almost the whole root, from the root trichoblast epidermis to vasculature, and from root tip to elongation and maturation zones [30], while depleted in the shoot [30] in control plants, as also found in the CK plants in the present study, that mRNAs were enriched for antioxidative defense system in the root comparing with the shoot. Unsurprisingly, mRNAs were also enriched for proliferation in the root. Taken together, the discrepancy of enriched gene function categories can be reasonably explained by the function differences between the two tissues.

Then transcriptomes of ‘Zhongshansa’ roots and shoots after 1 h of root waterlogging were compared with those under normal conditions. In total, there were 2310 (2.1%) DEGs found. Among them, 2090 DEGs were found in roots and 394 in shoots, indicating that the impact of soil waterlogging stress on ‘Zhongshansa’ transcripts was mainly in the stressed tissue, but that the aerial parts were also affected, as also shown in cotton [22] and Arabidopsis [32]. There were 174 DEGs shared by the two tissues, while the majority were tissue-specific. Consistent with this is that Marc et al. found the ability to tolerate hypoxic stress in roots and shoots could be genetically separate [7], and the anaerobic induction of most known ANPs were root specific in maize and Arabidopsis. It is not surprising that many of the tissue-specific DEGs were caused by the existence of tissue-specific cell populations, like photosynthetic cells in leaves. However, there are also some cell populations that have similar functions in both tissues, such asphloem cells, and their transcript changes under stress may contribute to the shared DEGs.

Effects on antioxidative defense system

Mustroph et al. compared transcriptomic adjustments to low-oxygen stress in 21 organisms across four kingdoms (Plantae, Animalia, Fungi and Bacteria) and found that the induction of enzymes that ameliorate ROS was a universal stress response, found in the majority of the evaluated species and especially in all plants [33]. When plants suffer from partial submergence, oxygen concentration in the root zone falls. With molecular oxygen being reduced to toxic reactive oxygen species (ROS) such as hydrogen peroxide, hydroxyl radicals, singlet oxygen and superoxide radicals [34], the balance between the production and quenching of the ROS in plants will be disrupted, which is critical to cell survival during flooding stress [35],[36]. To prevent the formation of ROS under stress, plants have evolved a complex antioxidative defense system: low molecular mass antioxidants (ascorbic acid, glutathione and tocopherols), enzymes regenerating the reduced forms of antioxidants, and ROS interacting enzymes such as superoxide dismutase, peroxidases and catalases [37]. Many antioxidant enzymes have been proven to be critical for many plants’ survival under different levels of waterlogging, e.g. tomato [34], eggplant [34], poplar [38], winter wheat [39], mungbean [40] and citrus [41]. The antioxidative defense system was greatly activated in CT root of ‘Zhongshansa’. Consistent with numerous studies that have shown a correlation between the ability to ameliorate ROS and survival under different levels of waterlogging, the high induction of ROS network proteins in waterlogged ‘Zhongshansa’ showed that strong detoxification was critical for its survival.

Effects on energy-consuming biosynthetic processes

Waterlogging led to a great repression in biogenesis of ribosomes, organelles and many other biosynthetic activities in ‘Zhongshansa’ roots. Notably, function categories related to cell proliferation were among the most enriched in CK root, while they were also dramatically depleted in treated root, which indicated a large scale of energy saving under hypoxic conditions. Energy-consuming biosynthesis processes of cellulose and cell wall were also greatly inhibited in CT shoot. Under waterlogging conditions, the mitochondrial respiration was inhibited and energy yield of alcoholic fermentation was significantly lower compared with respiration, which causes an energy crisis in anaerobic root [11]. The biological significance of a widespread inhibition of energy-consuming biosynthetic processes under waterlogging stress may be because it allows a concomitant reduction of ATP consumption [13]. Mustroph’s research showed that the restrictionof ATP-consuming processes like biogenesis of ribosomes, organelles and cell walls is an evolutionarily conserved coping mechanism across prokaryotes and eukaryotes [42]. The large-scale decline of mRNAs associated with biosynthetic processes in both tissues indicated that waterlogging of roots induces systemic inhibition of ATP-consuming processes.

Effects on carbon metabolism and amino acid metabolism

Since plants lack a circulatory system to mobilize oxygen produced by photosynthesis to heterotrophic roots [33], under waterlogging conditions the oxygen-dependent mitochondrial respiration in the root is greatly limited. A comparative analysis between plant species of transcriptional responses to hypoxia found contrasting expression profiles between the tolerant and susceptible species for genes encoding components of the mitochondrial electron transport chain, with genes mainly up-regulated in Arabidopsis, but down-regulated in poplar or rice [43]. In CT ‘Zhongshansa’ root, the majority of genes involved in the mitochondrial electron transport were down-regulated. Whether the mitochondrial electron transport chain transcript changes are related to plant waterlogging tolerance requires further demonstration.

As expected and verified by qRT-PCR, many genes including well-known hypoxic genes associated with glycolysis and fermentation (ADH, PDC and LDH were induced by waterlogging, indicated that the glycolysis and fermentation pathway was activated to maintain ATP production under the stress. As a result, the demand for carbohydrates increased, and significantly increased carbohydrate metabolism in treated roots. The acceleration of carbohydrate metabolism is conversed functionally among plants under hypoxic conditions, and has been proved to be critical for plants’ survival [33],[43]. Notably, in the shoot, two DEGs involved in starch and sucrose metabolism, annotated as encoding TPS, were up-regulated under the stress. The comparisons of early transcriptomes of poplar, Arabidopsis and cotton responses to waterlogging found that hypoxia triggers the overexpression of TPS in all three species [44]. TPS catalyzes the first step of trehalose synthesis, which is important in plant response to abiotic stresses [45]. TPS has been shown to regulate sugar metabolism in plants [46],[47], so the up-regulation of TPS in the shoot indicated the acceleration of sugar metabolism in the ‘Zhongshansa’ shoot. Many researchers have considered that shoots would transport carbohydrate to the root to supply more carbohydrates to hypoxic tissues, due to the higher demand for carbohydrates in glycolysis. In hypoxia-treated poplar, increased phloem transport of sucrose from leaves to roots was found [11]-research on Arabidopsis [32] and cotton [5] reached the same conclusion. So, the stimulation of starch and sucrose metabolismin ‘Zhongshansa’ shoot may also be involved in the systemic communications between anaerobic parts and aerial parts to survive soil waterlogging.

The comparative analysis of early transcriptome responses to low-oxygen environments in Arabidopsis, cotton and poplar found that amino acid metabolism changes were common in these three dicotyledonous species, although there was almost no overlap between their particular responses [44]. Waterlogging also led to rapid changes in the levels of amino acids in ‘Zhongshansa’ roots. In the CT, transcriptional down-regulation of genesinvolved in serine degradation was found. However, large numbers of genesinvolved in aspartic acid degradation were up-regulated. As a result, a rapid increase in serine and decrease in aspartic acid maybe found in the root. The same dynamic changes were found in the metabolite profiling of gray poplar root during hypoxia [11]. Kreuzwieser et al. proposed that hypoxia led to the inhibition of the TCA cycle and activation of glycolysis and fermentation pathways, resulting in an accumulation of amino acids closely derived from intermediates of glycolysis (e.g. serine) and a decrease of TCA cycle intermediate-derived amino acids (e.g. aspartic acid) [11].

Effects on non-symbiotic hemoglobins

Recent research by Narsai et al. on comparative analysis between plant species of transcriptional and metabolic responses to hypoxia paid special attention to the possible relationship between hemoglobin expression and plant tolerance to low-oxygen conditions [43],[48]. In plants, this protein family includes the symbiotic and non-symbiotic hemoglobins, the former are only expressed in nodules of legumes and some other species, and so the non-symbiotic hemoglobins are more commonly discussed in most plants. Narsai et al. found that transcript abundance of class-1 non-symbiotic hemoglobins rapidly increased under hypoxia in intolerant Arabidopsis, but were down-regulated or unchanged in tolerant rice and poplar; genes encoding class-2 and class-3 hemoglobins also showed similar but less extreme trends [43],[48]. However, an analysis of adaptive responses of two oak species to flooding stress suggested an inverse relationship between class-1 non-symbiotic hemoglobins gene expression and flooding tolerance [49]. Moreover, root transcript profiling analysis showed that submergence stress up-regulated hemoglobin in two flooding tolerant Rorippa species [50]. And expression of the gene encoding hemoglobin in cucumber susceptible to flooding stress decreased under waterlogging [20]. Parent et al. proposed that the interaction between non-symbiotic hemoglobins and nitric oxide (NO) was an alternative to the fermentation pathway under hypoxia, in which non-symbiotic hemoglobins acted as an NO dioxygenase to convert NO to nitrate [49],[51]. This pathway not only eliminated the toxic NO in the cell, but also helped maintain ATP synthesis. There were four unigenes annotated as encoding non-symbiotic hemoglobin: comp56074_c0, comp68739_c1, comp1346274_c0 and comp13327_c0. Among them, comp56074_c0 was annotated as encoding class 1 non-symbiotic hemoglobin and was 4.29-fold increased in CT root. According to the above hypothesis, this may facilitate ‘Zhongshansa’ waterlogging tolerance. Additionally, mRNAs related to heme (i.e. non-symbiotic hemoglobin) binding were highly up-regulated in the CT root. Thus, more attention needs to be paid to the relationship between non-symbiotic hemoglobins with ‘Zhongshansa’ waterlogging tolerance.

Conclusion

The transcript comparision of ‘Zhongshansa 406’ under normoxic and hypoxic condition using RNA-seq helped to explain the molecular basis of the early response of the remarkably waterlogging-tolerant ‘Zhongshansa’. Transcripts involved in the ROS network and carbon and nitrogen metabolism were greatly changed. In ‘Zhongshansa’, the former was mediated by the induction of active antioxidative defense system. To produce ATP, the glycolysis and fermentation pathway was stimulated and, as a result, sugars were supplied from the whole plant. Additionally, a series of ATP-consuming biosynthetic processes were dramatically repressed in shoots and roots. Amino acid metabolism was greatly changed, and a nitrate production pathway may also be induced to help maintain ATP. Thus transcript patterns reveal that ROS detoxification and energy maintenance were the primary coping mechanisms that ‘Zhongshansa’ adopted to survive oxygen deficiency, which may be responsible for its remarkable waterlogging tolerance.