- Email: [email protected]
Comprehensive analysis of differentially expressed genes reveals the molecular response to elevated CO2 levels in two sea buckthorn cultivars
Gene 660 (2018) 120–127
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Research paper
Comprehensive analysis of differentially expressed genes reveals the molecular response to elevated CO2 levels in two sea buckthorn cultivars ⁎
Guoyun Zhanga, Tong Zhanga, Juanjuan Liua, Jianguo Zhanga,b, , Caiyun Hea,
T
⁎
a State Key Laboratory of Tree Genetics and Breeding & Key Laboratory of Tree Breeding and Cultivation, State Forestry Administration, Research Institute of Forestry, Chinese Academy of Forestry, Beijing, China b Collaborative Innovation Center of Sustainable Forestry in Southern China, Nanjing Forestry University, Nanjing, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Sea buckthorn CO2 concentration Transcriptome Differentially expressed genes Molecular mechanism
Atmospheric carbon dioxide (CO2) concentration increases every year. It is critical to understand the elevated CO2 response molecular mechanisms of plants using genomic techniques. Hippophae rhamnoides L. is a high stress resistance plant species widely distributed in Europe and Asia. However, the molecular mechanism of elevated CO2 response in H. rhamnoides has been limited. In this study, transcriptomic analysis of two sea buckthorn cultivars under different CO2 concentrations was performed, based on the next-generation illumina sequencing platform and de novo assembly. We identified 4740 differentially expressed genes in sea buckthorn response to elevated CO2 concentrations. According to the gene ontology (GO) results, photosystem I, photosynthesis and chloroplast thylakoid membrane were the main enriched terms in ‘xiangyang’ sea buckthorn. In ‘zhongguo’ sea buckthorn, photosynthesis was also the main significantly enriched term. However, the number of photosynthesis related differentially expressed genes were different between two sea buckthorn cultivars. Our GO and pathway analyses indicated that the expression levels of the transcription factors WRKY, MYB and NAC were significantly different between the two sea buckthorn cultivars. This study provides a reliable transcriptome sequence resource and is a valuable resource for genetic and genomic researches for plants under high CO2 concentration in the future.
1. Introduction Sea buckthorn (Hippophae rhamnoides L.), a thorny, nitrogen-fixing, deciduous shrub, has been widely cultivated for its nutritional and medicinal properties (Suryakumar and Gupta, 2011; Teleszko et al., 2015; Li et al., 2017). Sea buckthorn is environmentally important in coastal areas, as it fixes sand in place, preventing wind erosion (Ledwood and Shimwell, 1971; Singh, 2005). Transcriptomic studies have shown that the sea buckthorn tolerates cold and freeze stress well (Chaudhary and Sharma, 2015). Owing to its tolerance of a wide range of temperatures, the native habitat of sea buckthorn extends across Europe and Asia, including China, India, Mongolia, Russia, Sweden, Finland, and Norway (Lian and Chen, 1999). Recent studies of this plant have explored such diverse topics as the bioactive compounds it produces, and the mechanisms causing fruit ripening (Ghangal et al., 2012; Teleszko et al., 2015; Li et al., 2017). However, the molecular response mechanisms of sea buckthorn to elevated atmospheric carbon dioxide remain unknown.
Atmospheric CO2 considered as the primary carbon source for plant photosynthesis and a major driver of climate change. Recently, studies of CO2 become an important component of plant researches. Compared with the lowest concentrations, the average global atmospheric concentration of CO2 was 394 μmol·mol−1 in 2012 (Meinshausen et al., 2011). This concentration are predicted to surpass 550 μmol·mol−1 by 2050 and reach 700 μmol·mol−1 by the end of 2100 (Leakey et al., 2009b). A large number of studies reveals that elevated CO2 significantly affect the productivity and fitness of plant (Idso and Kimball, 1997; Ainsworth and Long, 2005; Springer and Ward, 2007; Becklin et al., 2017). For different plant species, the capacity to elevated atmospheric CO2 relies on a series of physiological and biochemical processes. Hence, understanding the molecular mechanisms of tree growth can promote us to find certain key candidate genes, transcriptional factor and metabolism pathway. In past few decades, CO2 studies found that several physiological characters related to plant growth were affected by elevated atmospheric CO2 (Idso and Kimball, 1997; Ehlers et al., 2015; Becklin et al., 2017). Furthermore, a small number of
Abbreviations: DEGs, differentially expressed genes; RPKM, reads per kilobase per million reads; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; ZG, zhongguo; XY, xiangyang; PPI, protein–protein interaction; RSEM, RNA-seq by Expectation Maximization; STRING, Search Tool for the Retrieval of Interacting Genes ⁎ Corresponding authors. E-mail addresses: [email protected] (G. Zhang), [email protected] (J. Zhang), [email protected] (C. He). https://doi.org/10.1016/j.gene.2018.03.057 Received 24 December 2017; Received in revised form 5 March 2018; Accepted 16 March 2018 Available online 22 March 2018 0378-1119/ © 2018 Elsevier B.V. All rights reserved.
Gene 660 (2018) 120–127
G. Zhang et al.
determination and RNA extraction. Total RNA from each sample was isolated using the RNeasy Plus Universal Mini Kit (Qiagen, Germany) according to the manufacturer's instructions. The quality and quantity of total RNA was checked with an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Then, equal amounts of purification RNA from each separate tissue (leaf, stem, and root) of the plants under T0 treatment, and leaves of plants exposed to all three treatments (T0, T1, T2) were pooled to construct the cDNA libraries using the mRNA-seq sample preparation kit (Illumina, CA). Then, following the manufacture's recommendations, first-strand cDNA was synthesized from the generated short fragments using random hexamer primers. Then, second-strand cDNA synthesis was performed in a reaction mix. These cDNA fragments were purified using a QiaQuick PCR Purification Kit (Qiagen, Germany) and ligated to sequencing adapters. Suitable fragments were selected to create the final cDNA library by PCR amplification. The resulting 12 cDNA libraries were individually loaded into an illumina flow cell and sequenced on Illumina Hiseq2000 sequencing platform.
studies have provided insight into plant community dynamics (Ainsworth and Long, 2005; Leakey et al., 2009b; Norby et al., 2016) and large-scale transcriptomic studies have indicated that the expression levels of plant genes related to both photosynthesis (Kumar et al., 2017) and leaf development (Leakey et al., 2009a) are sensitive to elevated CO2 levels. For example, transcriptomic analyses found that a series of photosynthesis regulatory genes were apparently affected by elevated CO2 in Jatropha curcas L. (Kumar et al., 2017). Previous studies on elevated CO2 related molecular responses in plants concentrated on model species and few tree species and these studies mainly reveal secondary metabolism related genes during the plants development (Idso and Kimball, 1997; Tallis et al., 2010; Niu et al., 2016). Nowadays, with the high-speed development and widespread of next generation sequencing technologies, de novo RNA-sequencing and high throughput sequencing have promoted discovery and analysis of new genes in tree species. These approaches also provided insight for better understanding of transcriptional patterns and gene expression in plant development in response to diverse environments. Due to lack of genomic resources, few transcriptome studies for different tissues or development stages were performed for sea buckthorn. However, there has been little research on gene differential expression for sea buckthorn under elevated CO2 concentration. In this study, we exposed two sea buckthorn cultivars to three different levels of atmospheric CO2, including the current ambient concentration and two predicted future atmospheric CO2 concentrations (550 and 700 μmol·mol−1). Then, we used RNA-seq techniques to identify the DEGs among two cultivars at the three levels of CO2 exposure. A series of bioinformatics methods were used to determine the crucial differentially expressed genes affected by elevated CO2 concentration and related metabolic pathways in sea buckthorn. This study first provided temporal expression pattern of elevated CO2 response genes and comprehensive transcriptome sets for sea buckthorn.
2.3. De novo assembly and functional annotation We filtered the raw data to generate clean reads by removing adaptor, unknown nucleotides and low-quality reads. All subsequent analyses were performed on these high-quality clean reads. To obtain more comprehensive transcripts, we mixed the clean reads from different tissues. The mixed clean short reads were assembled using Trinity (Grabherr et al., 2011). Then, all assembled transcriptome sequences were compared with various sequence databases, including the NCBI non-redundant database (nr), NT, Swiss-Prot, and KOG. Based on these annotation, we annotated the assembled unigenes with the GO database using BLAST2GO with an E-value ≤ 1.0E−5 (Conesa et al., 2005). The KEGG pathway annotations were performed by KOBAS to the KEGG database. 2.4. Identification of DEGs
2. Material and methods
We generated six mRNA libraries representing the two sea buckthorn cultivars exposed to the three different CO2 concentrations. The raw sequencing data were produced using Illumina sequencing. After the raw data were generated and the data processing was completed, the high-quality reads were mapped to assembled transcriptome reference sequences using Bowtie (Wang et al., 2011). We used RNA-seq by Expectation Maximization (RSEM) to count the number of mapped reads for each gene (Li and Dewey, 2011). The expression level of each unigene was calculated as reads per kilobase per million mapped reads (RPKM) (Ward et al., 2012). After calculating gene expression, we standardized the RPKM data using the trimmed mean of m-values (TMM) approach. We used DEGseq package to identify DEGs (Wang et al., 2010). The resulting p-values were adjusted using the BenjaminiHochberg method of controlling the false discovery rate (FDR) (Shen et al., 2013). We defined DEGs in the same plant exposed to different levels of CO2 as significantly affected if the adjusted p-value was < 0.001, the FDR was < 0.01, and the |log2ratio| was > 1. Then, we performed a saturation analysis to test whether the number of detected genes increased proportional to the number of sequences. We used Pearson's correlation analysis to test for correlations between the parallel experiments (i.e. T0 vs T1, T0 vs T2, and T1 vs T2). A high correlation between parallel experiments would indicate the reproducibility of our results and show that our methods were stable. We further anatomized selected DEGs using integrated bioinformatics methods.
2.1. Plant material and experimental treatments On 11 March 2010, triennial “zhongguo” (ZG) and “xiangyang” (XY) sea buckthorn seedlings were planted in individual plastic pots (20 cm × 26 cm × 34 cm) containing a 5:3:2 mixture of clay soil, sand and peat. Eighteen randomly selected pots were moved into three controlled-environment chambers (AGC-2; Zhejiang University Electrical Equipment Factory, Hangzhou, China). Each chamber measured 3.5 m × 2.2 m × 3.2 m (length × width × height), with an average daytime active photosynthetic radiation of 800 ± 20 μmol·m−2 ·s−1, and a relative humidity of 65 ± 5%. Then, seedlings were exposed to three different CO2 concentrations for three months (Table 1), starting 25 June 2010. The CO2 concentrations in the chambers were continuously monitored by automatic CO2 detection systems. The pots were rotated inside the chambers weekly to minimize the effects of microclimatic variation within the chambers. 2.2. RNA isolation and sequencing The leaves, stems, and roots of all plants collected for biomass Table 1 The treatment of Hippophae rhamnoides sample.
T0 T1 T2
CO2 concentration(μmol/ mol)
Day temperature(°C)
Night temperature(°C)
385 ± 20 550 ± 20 720 ± 20
25 ± 1 28 ± 1 31 ± 1
20 ± 1 23 ± 1 26 ± 1
2.5. Differential expression analysis Based on the annotations of the reference transcriptome sequences, we got annotation of selected DEGs. We used the cluster algorithm of 121
Gene 660 (2018) 120–127
G. Zhang et al.
gene expression dynamics in java to identify unique DEGs expression profiles. Each profile group contained several DEGs with similar expression patterns. All subsequent analyses were performed on these profile groups. Then DEGs profile groups were used for further GO enrichment analysis with BLAST2GO. We performed pathway analysis on each profile group with KOBAS to determine the differences in metabolic processes among tissues. A p-value < 0.05 was set as thresholds for significant enrichment terms (Wu et al., 2006).
Table 2 Summary of sequencing output statistics of ‘xiangyang’ and ‘zhongguo’ sea buckthorn.
2.6. Constructing networks for significant DEGs We constructed protein-protein interaction (PPI) networks to better understand significant DEGs biological functions. We used Search Tool for the Retrieval of Interacting Genes (STRING) and cytoscape software to construct PPI networks between selected significant DEGs, based on input sequence, structure, and other data (Szklarczyk et al., 2011). Each PPI network was constructed with reference to database of known functions. Some network nodes, which had many connections, were considered as represent critical nodes (Yu et al., 2004). We statistically analyzed the distribution of the critical nodes in the constructed PPI networks to identify the center of the network.
Sample
Clean read
Clean nucleotides
GC (%)
Q30 (%)
XYL XYR XYS ZGL ZGR ZGS XY_T0 XY_T1 XY_T2 ZG_T0 ZG_T1 ZG_T2
23460504 24581860 36011500 25281576 26278119 31740290 15228890 15087552 11671999 13048249 11063727 13927626
2.35G 2.46G 3.6G 2.53G 2.63G 3.17G 1.52G 1.51G 1.17G 1.3G 1.11G 1.39G
42.22 42.37 42.46 42.16 41.78 41.98 41.82 42.1 42.37 41.79 42.04 42.09
91.21 90.92 90 89.76 91.11 89.86 92.65 92.56 91.96 90.98 89.87 93.52
functional KOG terms. The three largest categories of protein function across our unigenes were general functional prediction only, posttranslational modification, protein turnover, chaperon, and signal transduction (Supplementary File 2). We used Blast2GO to identify GO terms associated with the sea buckthorn transcripts. Sea buckthorn is not a model plant, so we were able to annotate only 27,185 (28.55%) unigenes with GO terms. The annotated unigenes were categorized as biological processes (71,721), cellular components (54,547) and molecular functions (33,043) (Supplementary File 3). To identify active biological pathways in the sea buckthorn, we mapped all assembled unigenes onto KEGG pathways with BLASTX. We assigned 4531 KO identifiers to 12,547 unigenes. We identified 199 KEGG pathways associated with at least 10 unigenes. The most common pathways were ribosome, protein processing in endoplasmic reticulum, and plant hormone signal transduction and the primary functional categories were organismal systems, metabolism, genetic information processing, environmental information processing, and cellular processes (Fig. 2). These results suggested that the sea buckthorn unigenes were associated with metabolism, environmental information, and genetic information.
2.7. Quantitative real-time PCR analysis validation of DEGs To validate the identified DEGs, eleven up- or down-regulated genes were randomly selected for RT-qPCR analyses using 18s-rRNA as the internal control gene. Total RNA was extracted from the sea buckthorn leaf tissue exposed to three levels of CO2 concentrations using a modified CTAB method. Based on the manufacturer's instructions, each qRTPCR reaction was performed on a LightCycler 480 II Real-time PCR Instrument (Roche, Swiss) using the QuantiFast® SYBR® Green PCR Kit (Qiagen, Germany). Primers were designed with Beacon Designer v7.0 software. The amplification was achieved by the following PCR program of first denaturation 95 °C for 10 min, then 40 cycles of denaturation at 95 °C for 10 s, and annealing and extension at 60 °C. The relative expression levels of the selected transcripts were normalized to 18s-rRNA and then calculated using the 2−ΔΔCt method. 3. Results
3.3. DEGs in the sea buckthorn leaf after exposure to elevated CO2 concentrations
3.1. Sequencing and assembly
We used RNA-seq data to identify the genes in the two sea buckthorn cultivars that were differentially expressed among the three CO2 concentrations. We obtained > 1.5 billion clean reads from each leaf sample. We next mapped the clean reads to the assembled transcriptome reference sequences. We considered DEGs significant if |log2Foldchange| > 1 and q-value < 0.05. A total of 1986 DEGs were filtered in ZG sea buckthorn under different treatments, including 701 up-regulated and 1285 down-regulated unigenes (Fig. 3). In XY sea buckthorn, there were 1442 up-regulated and 1312 down-regulated unigenes caused by elevated CO2 concentration. Then, we used k-means clustering to further analyze the DEGs. Unigenes with similar expression patterns were got together. In the ZG cultivar, the 1986 DEGs were recovered in eight subclusters (Supplementary File 4); in XY cultivar, the 2754 DEGs were recovered in 20 subclusters (Supplementary File 5).
To obtain a comprehensive and whole-plant transcriptome, we mixed the sequencing data generated from the different tissue samples collected from each cultivar (Table 2). After data cleaning and quality checking, 198,785 assembled clean reads with length ≥ 200 bp were remained. We then assembled the contigs to generate 95,204 unigenes with a mean length of 713 bp (Supplementary File 1). The transcriptome dataset has been submitted to the NCBI database (accession number SRP067785; Supplementary File 1). Our transcriptome was successfully sequenced, indicating that further bioinformatics analyses were possible. 3.2. Functional annotation and classification To validate and annotate the assembled unigenes, sequenced-based alignments were performed against four public databases using BLASTX with a significant e-value threshold of 1e−5. We then blasted all assembled sequences against the nr database, resulting in 26,047 significant BLAST hits (27.36% of all assembled sequences). We successfully mapped 23.43% of all unigenes against conserved domains in the Swiss-Prot database. We classified all 26,047 unigenes into 12,634 protein families from the Pfam database (Fig. 1). To predict and classify the possible functions of these families, we searched our assembled transcripts against the KOG database. We classify 11,731 unigenes into 25 KOG categories, comprising 2553
3.4. GO terms and pathway functional enrichment analysis We investigated the biological importance of the identified DEGs by performing a functional enrichment analysis with R tools. Only two subclusters of DEGs in the ZG cultivar exposed to high levels of CO2 were significantly enriched (enrichment of unigenes > 3 and p < 0.05). The DEGs in these two subclusters were categorized into 507 functional groups based on sequence homology. Among these, 307 122
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 1. Histogram of the 30 most abundant pfam domains of sea buckthorn unigenes.
significantly enriched. Our GO enrichment analysis indicated that the most significantly enriched GO terms for the DEGs in these seven subclusters were photosystem I, photosynthesis, chloroplast thylakoid membrane, and spermidine biosynthesis (Supplementary File 1). Several DEGs were involved in different aspects of photosynthesis, including photosynthesis, light harvesting, photosystem I, chloroplast
groups were related to biological processes and 113 groups were related to cellular components. Our GO enrichment analysis indicated that the most significantly enriched GO terms for the DEGs in the ZG cultivars exposed to high CO2 were photosynthesis, photosystem I, chitinase activity, and chitin catabolic process (Supplementary File 1). In the XY sea buckthorn cultivar, seven DEG subclusters were
Fig. 2. KEGG classification of the assembled unigenes. (A) cellular processes (B) environmental information processing (C) genetic information processing (D) metabolism (E) organismal systems. 123
Gene 660 (2018) 120–127
G. Zhang et al.
by the end of the century (Meinshausen et al., 2011). Past researches reveal that increased CO2 concentration directly affect plant growth in different species (Drake et al., 1997; Becklin et al., 2017; Drake et al., 2017). The most suitable intercellular CO2 concentration for the ZG sea buckthorn cultivar is 263 μmol·mol−1 (Drake et al., 1997). These changes in physiological characteristics prompted us to study the related molecular processes at the transcriptome level. So, we explored gene expression in the leaves of two sea buckthorn cultivars using RNAseq technology, with the aim of identifying the genes that responded to increased CO2 concentrations. In this study, a de novo transcriptome with sequences from 12 sea buckthorn samples were assembled. We totally identified 198,785 transcripts, including 95,204 unigenes. Using BLAST, we annotated 26,047 unigenes against nr database. Our functional bioinformatics analysis generated an abundance of unannotated unigenes, which could represent a gene pool specific to the sea buckthorn. This comprehensive analysis of gene expression in the sea buckthorn will be useful for future investigations of the genes and molecular mechanisms associated with exposure to elevated CO2 concentrations. Furthermore, we identified differentially expression genes between sea buckthorn cultivars exposed to different CO2 concentrations. The total number of genes in XYL1 vs XYL0 and XYL2 vs XYL1 were much higher than those in XYL2 vs XYL0 group. In XYL1 vs XYL0 and XYL2 vs XYL0 groups, most of the DEGs were up-regulated in sea buckthorn under elevated CO2 concentrations, suggesting that the complexity of the response to increased CO2 concentrations. Contrast to XY sea buckthorn, the total number of genes in ZGL3 vs ZGL1 was much higher than those in the other two groups. Interestingly, different from XY sea buckthorn, most of the DEGs were up-regulated in ZG sea buckthorn under elevated CO2 concentrations, which indicated that the different response mechanisms between two sea buckthorn cultivars. As many of the DEGs were different between the two cultivars, these cultivars probably employ different molecular mechanisms in response to elevated CO2 concentrations. In XY sea buckthorn, the most common significantly enriched GO terms were related to the spermidine biosynthetic process, photosynthesis, and catalytic activity. Spermidine is widely distributed in various plants and plays an important role in growth, age prevention, and resistance to abiotic stressors (Gonzalez et al., 2011; Alcázar and Tiburcio, 2014; Do et al., 2014; Kumar et al., 2015). For example, spermidine significantly enhanced the cold resistance of Citrus reticulata (Wang and Liu, 2015). Our results indicated that spermidine may play important role in elevated CO2 response process. Moreover, we identified 13 up-regulated DEGs associated with photosynthesis in the XY sea buckthorn. This up-regulation may imply an increased rate of photosynthesis. The most common enriched GO terms in the ZG cultivar were similar to these in the XY cultivar: photosynthesis, photosystem I, chitinase activity, and chitin catabolic process. DEGs encoding proteins involved in photosynthesis were identified, suggesting that elevated atmospheric CO2 causes a metabolic response in the sea buckthorn. Elevated CO2 concentrations influence photosynthesis in many plants (Long et al., 2004; Ramalho et al., 2013; Kumar et al., 2017). Increased levels of CO2 may modulate the expression of several genes related to photosynthesis and the Calvin-cycle (Fukayama et al., 2009; Kaiser et al., 2017; Noctor and Mhamdi, 2017). Our GO results identified several transcripts associated with photosynthetic photochemical reactions in the ZG sea buckthorn, which are the main components of the photosynthetic reaction center (Kumar et al., 2017). Seven of these genes were up-regulated and two were down-regulated in the ZG cultivars exposed to high CO2 concentrations. We also identified two genes encoding the chlorophyll a/b binding protein (Supplementary File 1). In the XY sea buckthorn cultivar, we identified 16 photosynthesis-related genes. Interestingly, all these genes were upregulated, including three related to the Photosystem I reaction center subunit and three related to the chlorophyll a/b binding protein (Supplementary File 1). Our data clearly indicated that the
Fig. 3. Distribution of DEGs between two treatments in XY and ZG sea buckthorn. ZGL0, ZGL1 and ZGL2: ZG sea buckthorn leaf under T0, T1 and T2 treatment. XYL0, XYL1 and XYL2: XY sea buckthorn leaf under T0, T1 and T2 treatment.
thylakoid membrane, and photosystem I reaction center. We used KOBAS to identify the significantly enriched DEG pathways. To identify these pathways, we used a hypergeometric test and adjusted the p-values using the Benjamini–Hochberg method. In the ZG sea buckthorn cultivar, 214 DEGs were significantly enriched in 24 pathways (corrected p < 0.05). The most common significant pathways were photosynthesis, phenylalanine metabolism, biosynthesis of secondary metabolites, and porphyrin and chlorophyll metabolism (Supplementary File 1). In XY sea buckthorn cultivar, 300 DEGs were significantly enriched in 56 pathways (corrected p < 0.05). The most common significant pathways were photosynthesis, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and flavonoid biosynthesis (Supplementary File 1). Lots of DEGs were recovered in pathways related to the biosynthesis of secondary metabolites, including phenylpropanoid, flavonoid, stilbenoid, diarylheptanoid, and gingerol. 3.5. The PPI network for significant DEGs Based on the protein interaction data from the STRING database (Szklarczyk et al., 2011), we generated a PPI network for the significant DEGs from each cultivar using Cytoscape. Statistical analysis of the two PPI networks identified the main hub genes (where edge > 30; Supplementary File 1). These hub genes likely play important roles in the DEG network when CO2 concentration increased. 3.6. Validation of DEGs among different tissues using qRT-PCR To confirm the reliability of the high-throughput sequencing data, we used qRT-PCR to measure gene expression. We selected 20 DEGs across the two sea buckthorn cultivars. Of these, 19 were validated and one was not detected. In XY sea buckthorn cultivar, our RNA-seq and qRT-PCR results for 12 genes were consistent and all 12 genes were upregulated. Eight DEGs were consistent between RNA-seq and qRTPCR (Fig. 4A). In ZG sea buckthorn cultivar, our RNA-seq and qRT-PCR results for three genes were consistent (Fig. 4B). The primers that we designed for each gene and for 18S-rRNA are summarized in Supplementary File 1. 4. Discussion Since 1959, atmospheric CO2 concentration has increased from 318 μmol·mol−1 to ~400 μmol·mol−1, and may reach 1000 μmol·mol−1 124
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 4. qRT-PCR validation of RNA-seq results. Two comparisons for each gene are shown in the figure. Upper panel: RNA-seq RPKM value; lower panel: fold changes by qRT-PCR measured in the same subjects used for RNA-seq analysis.
photochemical reactions related to photosynthesis were enhanced in both sea buckthorn cultivars exposed to high levels of CO2 (Fig. 5). We concluded that elevated CO2 significantly affected the expression of several photosynthesis-related genes in the sea buckthorn. Our results showed that the expression levels of several transcription factors (TFs) were significantly affected by elevated CO2
concentrations. In this study, 20 down-regulated and 27 up-regulated TFs were identified in sea buckthorn. Among these TFs, AP/ERF and zinc finger family have been shown to have crucial roles in plant responses to abiotic stresses. Also, genes encoding proteins in the same protein family are likely to have similar expression patterns. Interestingly, most of these TFs have higher expression in sea buckthorn 125
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 5. Model for the response to elevated CO2 concentrations in two sea buckthorn cultivars. Red shows up-regulated, green shows down-regulated and blue shows up- and down-regulated. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.gene.2018.03.057.
under T1 than T2, indicating that slightly elevated atmospheric CO2 increased the expression of TFs. The expression levels of down-regulated TFs in sea buckthorn under T1 were further decreased in sea buckthorn under T2. Consistent with previous studies (Shinde et al., 2015; Yousfi et al., 2016; Wei et al., 2017; Zhang et al., 2018), our results were showing that WRKY, MYB and NAC TFs have vital roles in the plant response to elevated CO2 (Supplementary File 1). Past researches reveal that abiotic stress can significantly affect transmembrane transport in plant (Haak et al., 2017). Transmembrane transport processes play important role in response to stress in plant by rebalancing cellular ion homeostasis. In this study, 21 transcripts involved in transmembrane transport functions were significantly upregulated in response to high CO2 concentrations (Supplementary File 1). For example, zinc transporter and the ABC transporter family were identified, suggesting that these transporters are key for the response to elevated CO2 in sea buckthorn (Fig. 5). Moreover, we found 37 downregulated transmembrane transport genes in the sea buckthorn. Among these genes, phosphate transporter is the most part, which are key for stress response in plant (Baek et al., 2017). These results indicated that phosphate transporter is involved in sea buckthorn response to increased CO2 concentrations. In total, our results indicated that elevated CO2 affected the expression of genes associated with transmembrane transport functions (Supplementary File 1).
Acknowledgement We would like to thank many other people for helping with the experiment. This study was supported by the Special Fund for Forest Scientific Research in the Public Welfare (201504103) and National Natural Science Foundation of China (31470616). Author contributions statement G. Y. Z, C.Y.H, and J.G.Z. performed the data analysis and drafted the manuscript. T.Z participated in the revision of manuscript. G. Y. Z, C.Y.H and J.J.L participated in the analysis of the data. C.Y.H and J.J.L performed the experiments. Competing financial interests The authors declare no competing financial interests. References Ainsworth, E.A., Long, S.P., 2005. What have we learned from 15 years of free-air CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2. New Phytol. 165, 351–371. Alcázar, R., Tiburcio, A.F., 2014. Plant polyamines in stress and development: an emerging area of research in plant sciences. Front. Plant Sci. 5, 319. Baek, D., Chun, H.J., Yun, D.J., Kim, M.C., 2017. Cross-talk between phosphate starvation and other environmental stress signaling pathways in plants. Mol. Cells 40, 697–705. Becklin, K.M., Walker 2nd, S.M., Way, D.A., Ward, J.K., 2017. CO2 studies remain key to understanding a future world. New Phytol. 214, 34–40. Chaudhary, S., Sharma, P.C., 2015. DeepSAGE Based Differential Gene Expression Analysis Under Cold and Freeze Stress in Seabuckthorn (Hippophae rhamnoides L.). Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. Do, P.T., Drechsel, O., Heyer, A.G., Hincha, D.K., Zuther, E., 2014. Changes in free polyamine levels, expression of polyamine biosynthesis genes, and performance of rice cultivars under salt stress: a comparison with responses to drought. Front. Plant Sci. 5. Drake, B.G., Gonzàlez-Meler, M.A., Long, S.P., 1997. More efficient plants: a consequence of rising atmospheric CO2? Annu. Rev. Plant Biol. 48, 609–639. Drake, B.L., Hanson, D.T., Lowrey, T.K., Sharp, Z.D., 2017. The carbon fertilization effect over a century of anthropogenic CO2 emissions: higher intracellular CO2 and more drought resistance among invasive and native grass species contrasts with increased water use efficiency for woody plants in the US Southwest. Glob. Chang. Biol. 23,
5. Conclusions The elevated atmospheric CO2 caused by global climate change has severe effect on plants. Here, to better understand the acclimation response of plants at the molecular level, we described changes in the transcriptome of the sea buckthorn leaf in response to various concentrations of CO2. Our RNA-seq analyses identified significant consistent changes in the expression of genes related to photosynthesis. Our GO and pathway analyses identified several significantly enriched photosynthesis related genes. The expression profiles of several genes and large-scale TFs were also altered in response to high levels of atmospheric CO2. Increased atmospheric CO2 caused many changes at the molecular level. Our results increase our understanding of the molecular mechanisms underlying the response of the sea buckthorn to elevated CO2. Our study provides a foundation for investigations of the molecular mechanisms that respond to elevated atmospheric CO2 concentrations in other plant species. 126
Gene 660 (2018) 120–127
G. Zhang et al.
Ellsworth, David S., Goll, Daniel S., Lapola, David M., Luus, Kristina A., Rob MacKenzie, A., Medlyn, Belinda E., Pavlick, Ryan, Rammig, Anja, Smith, Benjamin, Thomas, Rick, Thonicke, Kirsten, Walker, Anthony P., Yang, Xiaojuan, Zaehle, Sonke, 2016. Model–data synthesis for the next generation of forest free-air CO2 enrichment (FACE) experiments. New Phytol. 209, 17–28. Ramalho, J.C., Rodrigues, A.P., Semedo, J.N., Pais, I.P., Martins, L.D., Simões-Costa, M.C., Leitão, A.E., Fortunato, A.S., Batista-Santos, P., Palos, I.M., 2013. Sustained Photosynthetic Performance of Coffea spp. Under Long-term Enhanced [CO2]. Shen, Y., Zhang, Y., Chen, J., Lin, H., Zhao, M., Peng, H., Liu, L., Yuan, G., Zhang, S., Zhang, Z., 2013. Genome expression profile analysis reveals important transcripts in maize roots responding to the stress of heavy metal Pb. Physiol. Plant. 147, 270–282. Shinde, S., Behpouri, A., McElwain, J.C., Ng, C.K.-Y., 2015. Genome-wide transcriptomic analysis of the effects of sub-ambient atmospheric oxygen and elevated atmospheric carbon dioxide levels on gametophytes of the moss, Physcomitrella patens. J. Exp. Bot. erv197. Singh, V., 2005. In: Singh, V. (Ed.), Sea buckthorn (Hippophae L.) A Multipurpose Wonder Plant. vol. 2 Daya Publishing House. Springer, C.J., Ward, J.K., 2007. Flowering time and elevated atmospheric CO2. New Phytol. 176, 243–255. Suryakumar, G., Gupta, A., 2011. Medicinal and therapeutic potential of sea buckthorn (Hippophae rhamnoides L.). J. Ethnopharmacol. 138, 268–278. Szklarczyk, D., Franceschini, A., Kuhn, M., Simonovic, M., Roth, A., Minguez, P., Doerks, T., Stark, M., Muller, J., Bork, P., 2011. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 39, D561–D568. Tallis, M.J., Lin, Y., Rogers, A., Zhang, J., Street, N.R., Miglietta, F., Karnosky, D.F., De Angelis, P., Calfapietra, C., Taylor, G., 2010. The transcriptome of Populus in elevated CO reveals increased anthocyanin biosynthesis during delayed autumnal senescence. New Phytol. 186, 415–428. Teleszko, M., Wojdylo, A., Rudzinska, M., Oszmianski, J., Golis, T., 2015. Analysis of lipophilic and hydrophilic bioactive compounds content in sea buckthorn (Hippophae rhamnoides L.) berries. J. Agric. Food Chem. 63, 4120–4129. Wang, W., Liu, J.-H., 2015. Genome-wide identification and expression analysis of the polyamine oxidase gene family in sweet orange (Citrus sinensis). Gene 555, 421–429. Wang, L., Feng, Z., Wang, X., Wang, X., Zhang, X., 2010. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. Wang, W., Wei, Z., Lam, T.-W., Wang, J., 2011. Next generation sequencing has lower sequence coverage and poorer SNP-detection capability in the regulatory regions. Sci. Rep. 1. Ward, J.A., Ponnala, L., Weber, C.A., 2012. Strategies for transcriptome analysis in nonmodel plants. Am. J. Bot. 99, 267–276. Wei, Q., Luo, Q., Wang, R., Zhang, F., He, Y., Zhang, Y., Qiu, D., Li, K., Chang, J., Yang, G., He, G., 2017. A wheat R2R3-type MYB transcription factor TaODORANT1 positively regulates drought and salt stress responses in transgenic tobacco plants. Front. Plant Sci. 8, 1374. Wu, J., Mao, X., Cai, T., Luo, J., Wei, L., 2006. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 34, W720–W724. Yousfi, F.E., Makhloufi, E., Marande, W., Ghorbel, A.W., Bouzayen, M., Berges, H., 2016. Comparative analysis of WRKY genes potentially involved in salt stress responses in Triticum turgidum L. ssp. durum. Front. Plant Sci. 7, 2034. Yu, H., Greenbaum, D., Lu, H.X., Zhu, X., Gerstein, M., 2004. Genomic analysis of essentiality within protein networks. Trends Genet. 20, 227–231. Zhang, X., Yao, C., Fu, S., Xuan, H., Wen, S., Liu, C., Li, F., Liu, A., Bi, S., Zhang, S., Li, S., 2018. Stress2TF: a manually curated database of TF regulation in plant response to stress. Gene 638, 36–40.
782–792. Ehlers, I., Augusti, A., Betson, T.R., Nilsson, M.B., Marshall, J.D., Schleucher, J., 2015. Detecting long-term metabolic shifts using isotopomers: CO2-driven suppression of photorespiration in C3 plants over the 20th century. Proc. Natl. Acad. Sci. U. S. A. 112, 15585–15590. Fukayama, H., Fukuda, T., Masumoto, C., Taniguchi, Y., Sakai, H., Cheng, W., Hasegawa, T., Miyao, M., 2009. Rice plant response to long term CO2 enrichment: gene expression profiling. Plant Sci. 177, 203–210. Ghangal, R., Raghuvanshi, S., Sharma, P.C., 2012. Expressed sequence tag based identification and expression analysis of some cold inducible elements in seabuckthorn (Hippophae rhamnoides L.). Plant Physiol. Biochem. 51, 123–128. Gonzalez, M.E., Marco, F., Minguet, E.G., Sorli, P.C., Blázquez, M.A., Carbonell, J., Ruiz, O.A., Pieckenstain, F.L., 2011. Perturbation of spermine synthase gene expression and transcript profiling provide new insights on the role of the tetraamine spermine in Arabidopsis thaliana defense against Pseudomonas viridiflava. Plant Physiol. 110.171413. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. Haak, D.C., Fukao, T., Grene, R., Hua, Z., Ivanov, R., Perrella, G., Li, S., 2017. Multilevel regulation of abiotic stress responses in plants. Front. Plant Sci. 8, 1564. Idso, S., Kimball, B., 1997. Effects of long-term atmospheric CO2 enrichment on the growth and fruit production of sour orange trees. Glob. Chang. Biol. 3, 89–96. Kaiser, E., Zhou, D., Heuvelink, E., Harbinson, J., Morales, A., Marcelis, L.F.M., 2017. Elevated CO2 increases photosynthesis in fluctuating irradiance regardless of photosynthetic induction state. J. Exp. Bot. 68, 5629–5640. Kumar, M., Reddy, C., Ralph, P.J., 2015. Polyamines in morphogenesis and development: a promising research area in seaweeds. Front. Plant Sci. 6. Kumar, S., Sreeharsha, R.V., Mudalkar, S., Sarashetti, P.M., Reddy, A.R., 2017. Molecular insights into photosynthesis and carbohydrate metabolism in Jatropha curcas grown under elevated CO2 using transcriptome sequencing and assembly. Sci. Rep. 7, 11066. Leakey, A.D., Ainsworth, E.A., Bernacchi, C.J., Rogers, A., Long, S.P., Ort, D.R., 2009a. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876. Leakey, A.D., Ainsworth, E.A., Bernacchi, C.J., Rogers, A., Long, S.P., Ort, D.R., 2009b. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876. Ledwood, J., Shimwell, D., 1971. Growth-rates of Hippophaë rhamnoides L. Ann. Bot. 35, 1053–1058. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 12, 323. Li, Q.F., Zhang, Y.C., Chen, Y.Q., Yu, Y., 2017. Circular RNAs roll into the regulatory network of plants. Biochem. Biophys. Res. Commun. 488, 382–386. Lian, Y., Chen, X., 1999. The regular patterns of distribution on the natural components in plants of the genus Hippophae L. J. Northwest Norm. Univ. (Nat. Sci.) 36, 113–128. Long, S.P., Ainsworth, E.A., Rogers, A., Ort, D.R., 2004. Rising atmospheric carbon dioxide: plants FACE the future. Annu. Rev. Plant Biol. 55, 591–628. Meinshausen, M., Smith, S.J., Calvin, K., Daniel, J.S., Kainuma, M., Lamarque, J., Matsumoto, K., Montzka, S., Raper, S., Riahi, K., 2011. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. Clim. Chang. 109, 213–241. Niu, Y., Ahammed, G.J., Tang, C., Guo, L., Yu, J., 2016. Physiological and transcriptome responses to combinations of elevated CO2 and magnesium in Arabidopsis thaliana. PLoS One 11, e0149301. Noctor, G., Mhamdi, A., 2017. Climate change, CO2, and defense: the metabolic, redox, and signaling perspectives. Trends Plant Sci. 22, 857–870. Norby, Richard J., De Kauwe, M.G., Domingues, Tomas F., Duursma, Remko A.,
127
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Research paper
Comprehensive analysis of differentially expressed genes reveals the molecular response to elevated CO2 levels in two sea buckthorn cultivars ⁎
Guoyun Zhanga, Tong Zhanga, Juanjuan Liua, Jianguo Zhanga,b, , Caiyun Hea,
T
⁎
a State Key Laboratory of Tree Genetics and Breeding & Key Laboratory of Tree Breeding and Cultivation, State Forestry Administration, Research Institute of Forestry, Chinese Academy of Forestry, Beijing, China b Collaborative Innovation Center of Sustainable Forestry in Southern China, Nanjing Forestry University, Nanjing, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Sea buckthorn CO2 concentration Transcriptome Differentially expressed genes Molecular mechanism
Atmospheric carbon dioxide (CO2) concentration increases every year. It is critical to understand the elevated CO2 response molecular mechanisms of plants using genomic techniques. Hippophae rhamnoides L. is a high stress resistance plant species widely distributed in Europe and Asia. However, the molecular mechanism of elevated CO2 response in H. rhamnoides has been limited. In this study, transcriptomic analysis of two sea buckthorn cultivars under different CO2 concentrations was performed, based on the next-generation illumina sequencing platform and de novo assembly. We identified 4740 differentially expressed genes in sea buckthorn response to elevated CO2 concentrations. According to the gene ontology (GO) results, photosystem I, photosynthesis and chloroplast thylakoid membrane were the main enriched terms in ‘xiangyang’ sea buckthorn. In ‘zhongguo’ sea buckthorn, photosynthesis was also the main significantly enriched term. However, the number of photosynthesis related differentially expressed genes were different between two sea buckthorn cultivars. Our GO and pathway analyses indicated that the expression levels of the transcription factors WRKY, MYB and NAC were significantly different between the two sea buckthorn cultivars. This study provides a reliable transcriptome sequence resource and is a valuable resource for genetic and genomic researches for plants under high CO2 concentration in the future.
1. Introduction Sea buckthorn (Hippophae rhamnoides L.), a thorny, nitrogen-fixing, deciduous shrub, has been widely cultivated for its nutritional and medicinal properties (Suryakumar and Gupta, 2011; Teleszko et al., 2015; Li et al., 2017). Sea buckthorn is environmentally important in coastal areas, as it fixes sand in place, preventing wind erosion (Ledwood and Shimwell, 1971; Singh, 2005). Transcriptomic studies have shown that the sea buckthorn tolerates cold and freeze stress well (Chaudhary and Sharma, 2015). Owing to its tolerance of a wide range of temperatures, the native habitat of sea buckthorn extends across Europe and Asia, including China, India, Mongolia, Russia, Sweden, Finland, and Norway (Lian and Chen, 1999). Recent studies of this plant have explored such diverse topics as the bioactive compounds it produces, and the mechanisms causing fruit ripening (Ghangal et al., 2012; Teleszko et al., 2015; Li et al., 2017). However, the molecular response mechanisms of sea buckthorn to elevated atmospheric carbon dioxide remain unknown.
Atmospheric CO2 considered as the primary carbon source for plant photosynthesis and a major driver of climate change. Recently, studies of CO2 become an important component of plant researches. Compared with the lowest concentrations, the average global atmospheric concentration of CO2 was 394 μmol·mol−1 in 2012 (Meinshausen et al., 2011). This concentration are predicted to surpass 550 μmol·mol−1 by 2050 and reach 700 μmol·mol−1 by the end of 2100 (Leakey et al., 2009b). A large number of studies reveals that elevated CO2 significantly affect the productivity and fitness of plant (Idso and Kimball, 1997; Ainsworth and Long, 2005; Springer and Ward, 2007; Becklin et al., 2017). For different plant species, the capacity to elevated atmospheric CO2 relies on a series of physiological and biochemical processes. Hence, understanding the molecular mechanisms of tree growth can promote us to find certain key candidate genes, transcriptional factor and metabolism pathway. In past few decades, CO2 studies found that several physiological characters related to plant growth were affected by elevated atmospheric CO2 (Idso and Kimball, 1997; Ehlers et al., 2015; Becklin et al., 2017). Furthermore, a small number of
Abbreviations: DEGs, differentially expressed genes; RPKM, reads per kilobase per million reads; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; ZG, zhongguo; XY, xiangyang; PPI, protein–protein interaction; RSEM, RNA-seq by Expectation Maximization; STRING, Search Tool for the Retrieval of Interacting Genes ⁎ Corresponding authors. E-mail addresses: [email protected] (G. Zhang), [email protected] (J. Zhang), [email protected] (C. He). https://doi.org/10.1016/j.gene.2018.03.057 Received 24 December 2017; Received in revised form 5 March 2018; Accepted 16 March 2018 Available online 22 March 2018 0378-1119/ © 2018 Elsevier B.V. All rights reserved.
Gene 660 (2018) 120–127
G. Zhang et al.
determination and RNA extraction. Total RNA from each sample was isolated using the RNeasy Plus Universal Mini Kit (Qiagen, Germany) according to the manufacturer's instructions. The quality and quantity of total RNA was checked with an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Then, equal amounts of purification RNA from each separate tissue (leaf, stem, and root) of the plants under T0 treatment, and leaves of plants exposed to all three treatments (T0, T1, T2) were pooled to construct the cDNA libraries using the mRNA-seq sample preparation kit (Illumina, CA). Then, following the manufacture's recommendations, first-strand cDNA was synthesized from the generated short fragments using random hexamer primers. Then, second-strand cDNA synthesis was performed in a reaction mix. These cDNA fragments were purified using a QiaQuick PCR Purification Kit (Qiagen, Germany) and ligated to sequencing adapters. Suitable fragments were selected to create the final cDNA library by PCR amplification. The resulting 12 cDNA libraries were individually loaded into an illumina flow cell and sequenced on Illumina Hiseq2000 sequencing platform.
studies have provided insight into plant community dynamics (Ainsworth and Long, 2005; Leakey et al., 2009b; Norby et al., 2016) and large-scale transcriptomic studies have indicated that the expression levels of plant genes related to both photosynthesis (Kumar et al., 2017) and leaf development (Leakey et al., 2009a) are sensitive to elevated CO2 levels. For example, transcriptomic analyses found that a series of photosynthesis regulatory genes were apparently affected by elevated CO2 in Jatropha curcas L. (Kumar et al., 2017). Previous studies on elevated CO2 related molecular responses in plants concentrated on model species and few tree species and these studies mainly reveal secondary metabolism related genes during the plants development (Idso and Kimball, 1997; Tallis et al., 2010; Niu et al., 2016). Nowadays, with the high-speed development and widespread of next generation sequencing technologies, de novo RNA-sequencing and high throughput sequencing have promoted discovery and analysis of new genes in tree species. These approaches also provided insight for better understanding of transcriptional patterns and gene expression in plant development in response to diverse environments. Due to lack of genomic resources, few transcriptome studies for different tissues or development stages were performed for sea buckthorn. However, there has been little research on gene differential expression for sea buckthorn under elevated CO2 concentration. In this study, we exposed two sea buckthorn cultivars to three different levels of atmospheric CO2, including the current ambient concentration and two predicted future atmospheric CO2 concentrations (550 and 700 μmol·mol−1). Then, we used RNA-seq techniques to identify the DEGs among two cultivars at the three levels of CO2 exposure. A series of bioinformatics methods were used to determine the crucial differentially expressed genes affected by elevated CO2 concentration and related metabolic pathways in sea buckthorn. This study first provided temporal expression pattern of elevated CO2 response genes and comprehensive transcriptome sets for sea buckthorn.
2.3. De novo assembly and functional annotation We filtered the raw data to generate clean reads by removing adaptor, unknown nucleotides and low-quality reads. All subsequent analyses were performed on these high-quality clean reads. To obtain more comprehensive transcripts, we mixed the clean reads from different tissues. The mixed clean short reads were assembled using Trinity (Grabherr et al., 2011). Then, all assembled transcriptome sequences were compared with various sequence databases, including the NCBI non-redundant database (nr), NT, Swiss-Prot, and KOG. Based on these annotation, we annotated the assembled unigenes with the GO database using BLAST2GO with an E-value ≤ 1.0E−5 (Conesa et al., 2005). The KEGG pathway annotations were performed by KOBAS to the KEGG database. 2.4. Identification of DEGs
2. Material and methods
We generated six mRNA libraries representing the two sea buckthorn cultivars exposed to the three different CO2 concentrations. The raw sequencing data were produced using Illumina sequencing. After the raw data were generated and the data processing was completed, the high-quality reads were mapped to assembled transcriptome reference sequences using Bowtie (Wang et al., 2011). We used RNA-seq by Expectation Maximization (RSEM) to count the number of mapped reads for each gene (Li and Dewey, 2011). The expression level of each unigene was calculated as reads per kilobase per million mapped reads (RPKM) (Ward et al., 2012). After calculating gene expression, we standardized the RPKM data using the trimmed mean of m-values (TMM) approach. We used DEGseq package to identify DEGs (Wang et al., 2010). The resulting p-values were adjusted using the BenjaminiHochberg method of controlling the false discovery rate (FDR) (Shen et al., 2013). We defined DEGs in the same plant exposed to different levels of CO2 as significantly affected if the adjusted p-value was < 0.001, the FDR was < 0.01, and the |log2ratio| was > 1. Then, we performed a saturation analysis to test whether the number of detected genes increased proportional to the number of sequences. We used Pearson's correlation analysis to test for correlations between the parallel experiments (i.e. T0 vs T1, T0 vs T2, and T1 vs T2). A high correlation between parallel experiments would indicate the reproducibility of our results and show that our methods were stable. We further anatomized selected DEGs using integrated bioinformatics methods.
2.1. Plant material and experimental treatments On 11 March 2010, triennial “zhongguo” (ZG) and “xiangyang” (XY) sea buckthorn seedlings were planted in individual plastic pots (20 cm × 26 cm × 34 cm) containing a 5:3:2 mixture of clay soil, sand and peat. Eighteen randomly selected pots were moved into three controlled-environment chambers (AGC-2; Zhejiang University Electrical Equipment Factory, Hangzhou, China). Each chamber measured 3.5 m × 2.2 m × 3.2 m (length × width × height), with an average daytime active photosynthetic radiation of 800 ± 20 μmol·m−2 ·s−1, and a relative humidity of 65 ± 5%. Then, seedlings were exposed to three different CO2 concentrations for three months (Table 1), starting 25 June 2010. The CO2 concentrations in the chambers were continuously monitored by automatic CO2 detection systems. The pots were rotated inside the chambers weekly to minimize the effects of microclimatic variation within the chambers. 2.2. RNA isolation and sequencing The leaves, stems, and roots of all plants collected for biomass Table 1 The treatment of Hippophae rhamnoides sample.
T0 T1 T2
CO2 concentration(μmol/ mol)
Day temperature(°C)
Night temperature(°C)
385 ± 20 550 ± 20 720 ± 20
25 ± 1 28 ± 1 31 ± 1
20 ± 1 23 ± 1 26 ± 1
2.5. Differential expression analysis Based on the annotations of the reference transcriptome sequences, we got annotation of selected DEGs. We used the cluster algorithm of 121
Gene 660 (2018) 120–127
G. Zhang et al.
gene expression dynamics in java to identify unique DEGs expression profiles. Each profile group contained several DEGs with similar expression patterns. All subsequent analyses were performed on these profile groups. Then DEGs profile groups were used for further GO enrichment analysis with BLAST2GO. We performed pathway analysis on each profile group with KOBAS to determine the differences in metabolic processes among tissues. A p-value < 0.05 was set as thresholds for significant enrichment terms (Wu et al., 2006).
Table 2 Summary of sequencing output statistics of ‘xiangyang’ and ‘zhongguo’ sea buckthorn.
2.6. Constructing networks for significant DEGs We constructed protein-protein interaction (PPI) networks to better understand significant DEGs biological functions. We used Search Tool for the Retrieval of Interacting Genes (STRING) and cytoscape software to construct PPI networks between selected significant DEGs, based on input sequence, structure, and other data (Szklarczyk et al., 2011). Each PPI network was constructed with reference to database of known functions. Some network nodes, which had many connections, were considered as represent critical nodes (Yu et al., 2004). We statistically analyzed the distribution of the critical nodes in the constructed PPI networks to identify the center of the network.
Sample
Clean read
Clean nucleotides
GC (%)
Q30 (%)
XYL XYR XYS ZGL ZGR ZGS XY_T0 XY_T1 XY_T2 ZG_T0 ZG_T1 ZG_T2
23460504 24581860 36011500 25281576 26278119 31740290 15228890 15087552 11671999 13048249 11063727 13927626
2.35G 2.46G 3.6G 2.53G 2.63G 3.17G 1.52G 1.51G 1.17G 1.3G 1.11G 1.39G
42.22 42.37 42.46 42.16 41.78 41.98 41.82 42.1 42.37 41.79 42.04 42.09
91.21 90.92 90 89.76 91.11 89.86 92.65 92.56 91.96 90.98 89.87 93.52
functional KOG terms. The three largest categories of protein function across our unigenes were general functional prediction only, posttranslational modification, protein turnover, chaperon, and signal transduction (Supplementary File 2). We used Blast2GO to identify GO terms associated with the sea buckthorn transcripts. Sea buckthorn is not a model plant, so we were able to annotate only 27,185 (28.55%) unigenes with GO terms. The annotated unigenes were categorized as biological processes (71,721), cellular components (54,547) and molecular functions (33,043) (Supplementary File 3). To identify active biological pathways in the sea buckthorn, we mapped all assembled unigenes onto KEGG pathways with BLASTX. We assigned 4531 KO identifiers to 12,547 unigenes. We identified 199 KEGG pathways associated with at least 10 unigenes. The most common pathways were ribosome, protein processing in endoplasmic reticulum, and plant hormone signal transduction and the primary functional categories were organismal systems, metabolism, genetic information processing, environmental information processing, and cellular processes (Fig. 2). These results suggested that the sea buckthorn unigenes were associated with metabolism, environmental information, and genetic information.
2.7. Quantitative real-time PCR analysis validation of DEGs To validate the identified DEGs, eleven up- or down-regulated genes were randomly selected for RT-qPCR analyses using 18s-rRNA as the internal control gene. Total RNA was extracted from the sea buckthorn leaf tissue exposed to three levels of CO2 concentrations using a modified CTAB method. Based on the manufacturer's instructions, each qRTPCR reaction was performed on a LightCycler 480 II Real-time PCR Instrument (Roche, Swiss) using the QuantiFast® SYBR® Green PCR Kit (Qiagen, Germany). Primers were designed with Beacon Designer v7.0 software. The amplification was achieved by the following PCR program of first denaturation 95 °C for 10 min, then 40 cycles of denaturation at 95 °C for 10 s, and annealing and extension at 60 °C. The relative expression levels of the selected transcripts were normalized to 18s-rRNA and then calculated using the 2−ΔΔCt method. 3. Results
3.3. DEGs in the sea buckthorn leaf after exposure to elevated CO2 concentrations
3.1. Sequencing and assembly
We used RNA-seq data to identify the genes in the two sea buckthorn cultivars that were differentially expressed among the three CO2 concentrations. We obtained > 1.5 billion clean reads from each leaf sample. We next mapped the clean reads to the assembled transcriptome reference sequences. We considered DEGs significant if |log2Foldchange| > 1 and q-value < 0.05. A total of 1986 DEGs were filtered in ZG sea buckthorn under different treatments, including 701 up-regulated and 1285 down-regulated unigenes (Fig. 3). In XY sea buckthorn, there were 1442 up-regulated and 1312 down-regulated unigenes caused by elevated CO2 concentration. Then, we used k-means clustering to further analyze the DEGs. Unigenes with similar expression patterns were got together. In the ZG cultivar, the 1986 DEGs were recovered in eight subclusters (Supplementary File 4); in XY cultivar, the 2754 DEGs were recovered in 20 subclusters (Supplementary File 5).
To obtain a comprehensive and whole-plant transcriptome, we mixed the sequencing data generated from the different tissue samples collected from each cultivar (Table 2). After data cleaning and quality checking, 198,785 assembled clean reads with length ≥ 200 bp were remained. We then assembled the contigs to generate 95,204 unigenes with a mean length of 713 bp (Supplementary File 1). The transcriptome dataset has been submitted to the NCBI database (accession number SRP067785; Supplementary File 1). Our transcriptome was successfully sequenced, indicating that further bioinformatics analyses were possible. 3.2. Functional annotation and classification To validate and annotate the assembled unigenes, sequenced-based alignments were performed against four public databases using BLASTX with a significant e-value threshold of 1e−5. We then blasted all assembled sequences against the nr database, resulting in 26,047 significant BLAST hits (27.36% of all assembled sequences). We successfully mapped 23.43% of all unigenes against conserved domains in the Swiss-Prot database. We classified all 26,047 unigenes into 12,634 protein families from the Pfam database (Fig. 1). To predict and classify the possible functions of these families, we searched our assembled transcripts against the KOG database. We classify 11,731 unigenes into 25 KOG categories, comprising 2553
3.4. GO terms and pathway functional enrichment analysis We investigated the biological importance of the identified DEGs by performing a functional enrichment analysis with R tools. Only two subclusters of DEGs in the ZG cultivar exposed to high levels of CO2 were significantly enriched (enrichment of unigenes > 3 and p < 0.05). The DEGs in these two subclusters were categorized into 507 functional groups based on sequence homology. Among these, 307 122
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 1. Histogram of the 30 most abundant pfam domains of sea buckthorn unigenes.
significantly enriched. Our GO enrichment analysis indicated that the most significantly enriched GO terms for the DEGs in these seven subclusters were photosystem I, photosynthesis, chloroplast thylakoid membrane, and spermidine biosynthesis (Supplementary File 1). Several DEGs were involved in different aspects of photosynthesis, including photosynthesis, light harvesting, photosystem I, chloroplast
groups were related to biological processes and 113 groups were related to cellular components. Our GO enrichment analysis indicated that the most significantly enriched GO terms for the DEGs in the ZG cultivars exposed to high CO2 were photosynthesis, photosystem I, chitinase activity, and chitin catabolic process (Supplementary File 1). In the XY sea buckthorn cultivar, seven DEG subclusters were
Fig. 2. KEGG classification of the assembled unigenes. (A) cellular processes (B) environmental information processing (C) genetic information processing (D) metabolism (E) organismal systems. 123
Gene 660 (2018) 120–127
G. Zhang et al.
by the end of the century (Meinshausen et al., 2011). Past researches reveal that increased CO2 concentration directly affect plant growth in different species (Drake et al., 1997; Becklin et al., 2017; Drake et al., 2017). The most suitable intercellular CO2 concentration for the ZG sea buckthorn cultivar is 263 μmol·mol−1 (Drake et al., 1997). These changes in physiological characteristics prompted us to study the related molecular processes at the transcriptome level. So, we explored gene expression in the leaves of two sea buckthorn cultivars using RNAseq technology, with the aim of identifying the genes that responded to increased CO2 concentrations. In this study, a de novo transcriptome with sequences from 12 sea buckthorn samples were assembled. We totally identified 198,785 transcripts, including 95,204 unigenes. Using BLAST, we annotated 26,047 unigenes against nr database. Our functional bioinformatics analysis generated an abundance of unannotated unigenes, which could represent a gene pool specific to the sea buckthorn. This comprehensive analysis of gene expression in the sea buckthorn will be useful for future investigations of the genes and molecular mechanisms associated with exposure to elevated CO2 concentrations. Furthermore, we identified differentially expression genes between sea buckthorn cultivars exposed to different CO2 concentrations. The total number of genes in XYL1 vs XYL0 and XYL2 vs XYL1 were much higher than those in XYL2 vs XYL0 group. In XYL1 vs XYL0 and XYL2 vs XYL0 groups, most of the DEGs were up-regulated in sea buckthorn under elevated CO2 concentrations, suggesting that the complexity of the response to increased CO2 concentrations. Contrast to XY sea buckthorn, the total number of genes in ZGL3 vs ZGL1 was much higher than those in the other two groups. Interestingly, different from XY sea buckthorn, most of the DEGs were up-regulated in ZG sea buckthorn under elevated CO2 concentrations, which indicated that the different response mechanisms between two sea buckthorn cultivars. As many of the DEGs were different between the two cultivars, these cultivars probably employ different molecular mechanisms in response to elevated CO2 concentrations. In XY sea buckthorn, the most common significantly enriched GO terms were related to the spermidine biosynthetic process, photosynthesis, and catalytic activity. Spermidine is widely distributed in various plants and plays an important role in growth, age prevention, and resistance to abiotic stressors (Gonzalez et al., 2011; Alcázar and Tiburcio, 2014; Do et al., 2014; Kumar et al., 2015). For example, spermidine significantly enhanced the cold resistance of Citrus reticulata (Wang and Liu, 2015). Our results indicated that spermidine may play important role in elevated CO2 response process. Moreover, we identified 13 up-regulated DEGs associated with photosynthesis in the XY sea buckthorn. This up-regulation may imply an increased rate of photosynthesis. The most common enriched GO terms in the ZG cultivar were similar to these in the XY cultivar: photosynthesis, photosystem I, chitinase activity, and chitin catabolic process. DEGs encoding proteins involved in photosynthesis were identified, suggesting that elevated atmospheric CO2 causes a metabolic response in the sea buckthorn. Elevated CO2 concentrations influence photosynthesis in many plants (Long et al., 2004; Ramalho et al., 2013; Kumar et al., 2017). Increased levels of CO2 may modulate the expression of several genes related to photosynthesis and the Calvin-cycle (Fukayama et al., 2009; Kaiser et al., 2017; Noctor and Mhamdi, 2017). Our GO results identified several transcripts associated with photosynthetic photochemical reactions in the ZG sea buckthorn, which are the main components of the photosynthetic reaction center (Kumar et al., 2017). Seven of these genes were up-regulated and two were down-regulated in the ZG cultivars exposed to high CO2 concentrations. We also identified two genes encoding the chlorophyll a/b binding protein (Supplementary File 1). In the XY sea buckthorn cultivar, we identified 16 photosynthesis-related genes. Interestingly, all these genes were upregulated, including three related to the Photosystem I reaction center subunit and three related to the chlorophyll a/b binding protein (Supplementary File 1). Our data clearly indicated that the
Fig. 3. Distribution of DEGs between two treatments in XY and ZG sea buckthorn. ZGL0, ZGL1 and ZGL2: ZG sea buckthorn leaf under T0, T1 and T2 treatment. XYL0, XYL1 and XYL2: XY sea buckthorn leaf under T0, T1 and T2 treatment.
thylakoid membrane, and photosystem I reaction center. We used KOBAS to identify the significantly enriched DEG pathways. To identify these pathways, we used a hypergeometric test and adjusted the p-values using the Benjamini–Hochberg method. In the ZG sea buckthorn cultivar, 214 DEGs were significantly enriched in 24 pathways (corrected p < 0.05). The most common significant pathways were photosynthesis, phenylalanine metabolism, biosynthesis of secondary metabolites, and porphyrin and chlorophyll metabolism (Supplementary File 1). In XY sea buckthorn cultivar, 300 DEGs were significantly enriched in 56 pathways (corrected p < 0.05). The most common significant pathways were photosynthesis, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and flavonoid biosynthesis (Supplementary File 1). Lots of DEGs were recovered in pathways related to the biosynthesis of secondary metabolites, including phenylpropanoid, flavonoid, stilbenoid, diarylheptanoid, and gingerol. 3.5. The PPI network for significant DEGs Based on the protein interaction data from the STRING database (Szklarczyk et al., 2011), we generated a PPI network for the significant DEGs from each cultivar using Cytoscape. Statistical analysis of the two PPI networks identified the main hub genes (where edge > 30; Supplementary File 1). These hub genes likely play important roles in the DEG network when CO2 concentration increased. 3.6. Validation of DEGs among different tissues using qRT-PCR To confirm the reliability of the high-throughput sequencing data, we used qRT-PCR to measure gene expression. We selected 20 DEGs across the two sea buckthorn cultivars. Of these, 19 were validated and one was not detected. In XY sea buckthorn cultivar, our RNA-seq and qRT-PCR results for 12 genes were consistent and all 12 genes were upregulated. Eight DEGs were consistent between RNA-seq and qRTPCR (Fig. 4A). In ZG sea buckthorn cultivar, our RNA-seq and qRT-PCR results for three genes were consistent (Fig. 4B). The primers that we designed for each gene and for 18S-rRNA are summarized in Supplementary File 1. 4. Discussion Since 1959, atmospheric CO2 concentration has increased from 318 μmol·mol−1 to ~400 μmol·mol−1, and may reach 1000 μmol·mol−1 124
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 4. qRT-PCR validation of RNA-seq results. Two comparisons for each gene are shown in the figure. Upper panel: RNA-seq RPKM value; lower panel: fold changes by qRT-PCR measured in the same subjects used for RNA-seq analysis.
photochemical reactions related to photosynthesis were enhanced in both sea buckthorn cultivars exposed to high levels of CO2 (Fig. 5). We concluded that elevated CO2 significantly affected the expression of several photosynthesis-related genes in the sea buckthorn. Our results showed that the expression levels of several transcription factors (TFs) were significantly affected by elevated CO2
concentrations. In this study, 20 down-regulated and 27 up-regulated TFs were identified in sea buckthorn. Among these TFs, AP/ERF and zinc finger family have been shown to have crucial roles in plant responses to abiotic stresses. Also, genes encoding proteins in the same protein family are likely to have similar expression patterns. Interestingly, most of these TFs have higher expression in sea buckthorn 125
Gene 660 (2018) 120–127
G. Zhang et al.
Fig. 5. Model for the response to elevated CO2 concentrations in two sea buckthorn cultivars. Red shows up-regulated, green shows down-regulated and blue shows up- and down-regulated. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.gene.2018.03.057.
under T1 than T2, indicating that slightly elevated atmospheric CO2 increased the expression of TFs. The expression levels of down-regulated TFs in sea buckthorn under T1 were further decreased in sea buckthorn under T2. Consistent with previous studies (Shinde et al., 2015; Yousfi et al., 2016; Wei et al., 2017; Zhang et al., 2018), our results were showing that WRKY, MYB and NAC TFs have vital roles in the plant response to elevated CO2 (Supplementary File 1). Past researches reveal that abiotic stress can significantly affect transmembrane transport in plant (Haak et al., 2017). Transmembrane transport processes play important role in response to stress in plant by rebalancing cellular ion homeostasis. In this study, 21 transcripts involved in transmembrane transport functions were significantly upregulated in response to high CO2 concentrations (Supplementary File 1). For example, zinc transporter and the ABC transporter family were identified, suggesting that these transporters are key for the response to elevated CO2 in sea buckthorn (Fig. 5). Moreover, we found 37 downregulated transmembrane transport genes in the sea buckthorn. Among these genes, phosphate transporter is the most part, which are key for stress response in plant (Baek et al., 2017). These results indicated that phosphate transporter is involved in sea buckthorn response to increased CO2 concentrations. In total, our results indicated that elevated CO2 affected the expression of genes associated with transmembrane transport functions (Supplementary File 1).
Acknowledgement We would like to thank many other people for helping with the experiment. This study was supported by the Special Fund for Forest Scientific Research in the Public Welfare (201504103) and National Natural Science Foundation of China (31470616). Author contributions statement G. Y. Z, C.Y.H, and J.G.Z. performed the data analysis and drafted the manuscript. T.Z participated in the revision of manuscript. G. Y. Z, C.Y.H and J.J.L participated in the analysis of the data. C.Y.H and J.J.L performed the experiments. Competing financial interests The authors declare no competing financial interests. References Ainsworth, E.A., Long, S.P., 2005. What have we learned from 15 years of free-air CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2. New Phytol. 165, 351–371. Alcázar, R., Tiburcio, A.F., 2014. Plant polyamines in stress and development: an emerging area of research in plant sciences. Front. Plant Sci. 5, 319. Baek, D., Chun, H.J., Yun, D.J., Kim, M.C., 2017. Cross-talk between phosphate starvation and other environmental stress signaling pathways in plants. Mol. Cells 40, 697–705. Becklin, K.M., Walker 2nd, S.M., Way, D.A., Ward, J.K., 2017. CO2 studies remain key to understanding a future world. New Phytol. 214, 34–40. Chaudhary, S., Sharma, P.C., 2015. DeepSAGE Based Differential Gene Expression Analysis Under Cold and Freeze Stress in Seabuckthorn (Hippophae rhamnoides L.). Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. Do, P.T., Drechsel, O., Heyer, A.G., Hincha, D.K., Zuther, E., 2014. Changes in free polyamine levels, expression of polyamine biosynthesis genes, and performance of rice cultivars under salt stress: a comparison with responses to drought. Front. Plant Sci. 5. Drake, B.G., Gonzàlez-Meler, M.A., Long, S.P., 1997. More efficient plants: a consequence of rising atmospheric CO2? Annu. Rev. Plant Biol. 48, 609–639. Drake, B.L., Hanson, D.T., Lowrey, T.K., Sharp, Z.D., 2017. The carbon fertilization effect over a century of anthropogenic CO2 emissions: higher intracellular CO2 and more drought resistance among invasive and native grass species contrasts with increased water use efficiency for woody plants in the US Southwest. Glob. Chang. Biol. 23,
5. Conclusions The elevated atmospheric CO2 caused by global climate change has severe effect on plants. Here, to better understand the acclimation response of plants at the molecular level, we described changes in the transcriptome of the sea buckthorn leaf in response to various concentrations of CO2. Our RNA-seq analyses identified significant consistent changes in the expression of genes related to photosynthesis. Our GO and pathway analyses identified several significantly enriched photosynthesis related genes. The expression profiles of several genes and large-scale TFs were also altered in response to high levels of atmospheric CO2. Increased atmospheric CO2 caused many changes at the molecular level. Our results increase our understanding of the molecular mechanisms underlying the response of the sea buckthorn to elevated CO2. Our study provides a foundation for investigations of the molecular mechanisms that respond to elevated atmospheric CO2 concentrations in other plant species. 126
Gene 660 (2018) 120–127
G. Zhang et al.
Ellsworth, David S., Goll, Daniel S., Lapola, David M., Luus, Kristina A., Rob MacKenzie, A., Medlyn, Belinda E., Pavlick, Ryan, Rammig, Anja, Smith, Benjamin, Thomas, Rick, Thonicke, Kirsten, Walker, Anthony P., Yang, Xiaojuan, Zaehle, Sonke, 2016. Model–data synthesis for the next generation of forest free-air CO2 enrichment (FACE) experiments. New Phytol. 209, 17–28. Ramalho, J.C., Rodrigues, A.P., Semedo, J.N., Pais, I.P., Martins, L.D., Simões-Costa, M.C., Leitão, A.E., Fortunato, A.S., Batista-Santos, P., Palos, I.M., 2013. Sustained Photosynthetic Performance of Coffea spp. Under Long-term Enhanced [CO2]. Shen, Y., Zhang, Y., Chen, J., Lin, H., Zhao, M., Peng, H., Liu, L., Yuan, G., Zhang, S., Zhang, Z., 2013. Genome expression profile analysis reveals important transcripts in maize roots responding to the stress of heavy metal Pb. Physiol. Plant. 147, 270–282. Shinde, S., Behpouri, A., McElwain, J.C., Ng, C.K.-Y., 2015. Genome-wide transcriptomic analysis of the effects of sub-ambient atmospheric oxygen and elevated atmospheric carbon dioxide levels on gametophytes of the moss, Physcomitrella patens. J. Exp. Bot. erv197. Singh, V., 2005. In: Singh, V. (Ed.), Sea buckthorn (Hippophae L.) A Multipurpose Wonder Plant. vol. 2 Daya Publishing House. Springer, C.J., Ward, J.K., 2007. Flowering time and elevated atmospheric CO2. New Phytol. 176, 243–255. Suryakumar, G., Gupta, A., 2011. Medicinal and therapeutic potential of sea buckthorn (Hippophae rhamnoides L.). J. Ethnopharmacol. 138, 268–278. Szklarczyk, D., Franceschini, A., Kuhn, M., Simonovic, M., Roth, A., Minguez, P., Doerks, T., Stark, M., Muller, J., Bork, P., 2011. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 39, D561–D568. Tallis, M.J., Lin, Y., Rogers, A., Zhang, J., Street, N.R., Miglietta, F., Karnosky, D.F., De Angelis, P., Calfapietra, C., Taylor, G., 2010. The transcriptome of Populus in elevated CO reveals increased anthocyanin biosynthesis during delayed autumnal senescence. New Phytol. 186, 415–428. Teleszko, M., Wojdylo, A., Rudzinska, M., Oszmianski, J., Golis, T., 2015. Analysis of lipophilic and hydrophilic bioactive compounds content in sea buckthorn (Hippophae rhamnoides L.) berries. J. Agric. Food Chem. 63, 4120–4129. Wang, W., Liu, J.-H., 2015. Genome-wide identification and expression analysis of the polyamine oxidase gene family in sweet orange (Citrus sinensis). Gene 555, 421–429. Wang, L., Feng, Z., Wang, X., Wang, X., Zhang, X., 2010. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. Wang, W., Wei, Z., Lam, T.-W., Wang, J., 2011. Next generation sequencing has lower sequence coverage and poorer SNP-detection capability in the regulatory regions. Sci. Rep. 1. Ward, J.A., Ponnala, L., Weber, C.A., 2012. Strategies for transcriptome analysis in nonmodel plants. Am. J. Bot. 99, 267–276. Wei, Q., Luo, Q., Wang, R., Zhang, F., He, Y., Zhang, Y., Qiu, D., Li, K., Chang, J., Yang, G., He, G., 2017. A wheat R2R3-type MYB transcription factor TaODORANT1 positively regulates drought and salt stress responses in transgenic tobacco plants. Front. Plant Sci. 8, 1374. Wu, J., Mao, X., Cai, T., Luo, J., Wei, L., 2006. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 34, W720–W724. Yousfi, F.E., Makhloufi, E., Marande, W., Ghorbel, A.W., Bouzayen, M., Berges, H., 2016. Comparative analysis of WRKY genes potentially involved in salt stress responses in Triticum turgidum L. ssp. durum. Front. Plant Sci. 7, 2034. Yu, H., Greenbaum, D., Lu, H.X., Zhu, X., Gerstein, M., 2004. Genomic analysis of essentiality within protein networks. Trends Genet. 20, 227–231. Zhang, X., Yao, C., Fu, S., Xuan, H., Wen, S., Liu, C., Li, F., Liu, A., Bi, S., Zhang, S., Li, S., 2018. Stress2TF: a manually curated database of TF regulation in plant response to stress. Gene 638, 36–40.
782–792. Ehlers, I., Augusti, A., Betson, T.R., Nilsson, M.B., Marshall, J.D., Schleucher, J., 2015. Detecting long-term metabolic shifts using isotopomers: CO2-driven suppression of photorespiration in C3 plants over the 20th century. Proc. Natl. Acad. Sci. U. S. A. 112, 15585–15590. Fukayama, H., Fukuda, T., Masumoto, C., Taniguchi, Y., Sakai, H., Cheng, W., Hasegawa, T., Miyao, M., 2009. Rice plant response to long term CO2 enrichment: gene expression profiling. Plant Sci. 177, 203–210. Ghangal, R., Raghuvanshi, S., Sharma, P.C., 2012. Expressed sequence tag based identification and expression analysis of some cold inducible elements in seabuckthorn (Hippophae rhamnoides L.). Plant Physiol. Biochem. 51, 123–128. Gonzalez, M.E., Marco, F., Minguet, E.G., Sorli, P.C., Blázquez, M.A., Carbonell, J., Ruiz, O.A., Pieckenstain, F.L., 2011. Perturbation of spermine synthase gene expression and transcript profiling provide new insights on the role of the tetraamine spermine in Arabidopsis thaliana defense against Pseudomonas viridiflava. Plant Physiol. 110.171413. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. Haak, D.C., Fukao, T., Grene, R., Hua, Z., Ivanov, R., Perrella, G., Li, S., 2017. Multilevel regulation of abiotic stress responses in plants. Front. Plant Sci. 8, 1564. Idso, S., Kimball, B., 1997. Effects of long-term atmospheric CO2 enrichment on the growth and fruit production of sour orange trees. Glob. Chang. Biol. 3, 89–96. Kaiser, E., Zhou, D., Heuvelink, E., Harbinson, J., Morales, A., Marcelis, L.F.M., 2017. Elevated CO2 increases photosynthesis in fluctuating irradiance regardless of photosynthetic induction state. J. Exp. Bot. 68, 5629–5640. Kumar, M., Reddy, C., Ralph, P.J., 2015. Polyamines in morphogenesis and development: a promising research area in seaweeds. Front. Plant Sci. 6. Kumar, S., Sreeharsha, R.V., Mudalkar, S., Sarashetti, P.M., Reddy, A.R., 2017. Molecular insights into photosynthesis and carbohydrate metabolism in Jatropha curcas grown under elevated CO2 using transcriptome sequencing and assembly. Sci. Rep. 7, 11066. Leakey, A.D., Ainsworth, E.A., Bernacchi, C.J., Rogers, A., Long, S.P., Ort, D.R., 2009a. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876. Leakey, A.D., Ainsworth, E.A., Bernacchi, C.J., Rogers, A., Long, S.P., Ort, D.R., 2009b. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876. Ledwood, J., Shimwell, D., 1971. Growth-rates of Hippophaë rhamnoides L. Ann. Bot. 35, 1053–1058. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 12, 323. Li, Q.F., Zhang, Y.C., Chen, Y.Q., Yu, Y., 2017. Circular RNAs roll into the regulatory network of plants. Biochem. Biophys. Res. Commun. 488, 382–386. Lian, Y., Chen, X., 1999. The regular patterns of distribution on the natural components in plants of the genus Hippophae L. J. Northwest Norm. Univ. (Nat. Sci.) 36, 113–128. Long, S.P., Ainsworth, E.A., Rogers, A., Ort, D.R., 2004. Rising atmospheric carbon dioxide: plants FACE the future. Annu. Rev. Plant Biol. 55, 591–628. Meinshausen, M., Smith, S.J., Calvin, K., Daniel, J.S., Kainuma, M., Lamarque, J., Matsumoto, K., Montzka, S., Raper, S., Riahi, K., 2011. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. Clim. Chang. 109, 213–241. Niu, Y., Ahammed, G.J., Tang, C., Guo, L., Yu, J., 2016. Physiological and transcriptome responses to combinations of elevated CO2 and magnesium in Arabidopsis thaliana. PLoS One 11, e0149301. Noctor, G., Mhamdi, A., 2017. Climate change, CO2, and defense: the metabolic, redox, and signaling perspectives. Trends Plant Sci. 22, 857–870. Norby, Richard J., De Kauwe, M.G., Domingues, Tomas F., Duursma, Remko A.,
127