Transcriptome analysis identifies differentially expressed genes in maize leaf tissues in response to elevated atmospheric [CO2]

ABSTRACT The sustained increase in atmospheric [CO2] over the past two centuries has brought a series of global challenges in plant-environment interactions. However, genetic mechanisms in botanical adaption and feedback to environmental alteration remain elusive. Here we collected and analysed transcriptome sequencing data from leaf tissues of maize Zea mays conditioned under high [CO2] stress for 0, 7 or 14 d. A total of 1390 genes, either up- or down-regulated, differed significantly in these conditions. Gene Ontology (GO) enrichment terms and KEGG metabolism pathways included protein phosphorylation, protein ubiquitination, oxidation-reduction, plant organ development, cellular response to endogenous response and MAPK signaling. We also identified three densely connected gene clusters: TCP/PYE, WRKYs and MYC/JAZ from significant DEGs. In conclusion, we provided a genetic database for biologists in botany, agronomy and ecology to further investigate the molecular machinery by which C4-plants respond to elevated atmospheric [CO2] and accompanying global climate change.


Introduction
Carbon dioxide, an important trace component in the atmosphere, acts as an integral part of the global carbon cycle. Human activities, primarily the burning of fossil fuels and deforestation, have caused a remarkable and continuous increase in atmospheric [CO 2 ] over the past two centuries. According to NASA reports, global annual mean [CO 2 ] has increased from 280 ppm in the mid-18th century to 407 ppm in 2017 (Earth System Research Laboratory, NASA. Website: www.esrl.noaa.gov/gmd/ccgg/trends/). Consequently, the [CO 2 ] increase may result in intensified global warming, dubbed 'greenhouse effect'.
Plants consume solar energy, atmospheric CO 2 and water to produce carbohydrate in photosynthesis. The potential impacts of elevated [CO 2 ] in the air and accompanying climate change on land plants have attracted scientific attention since the 1990s. It is well known that elevated [CO 2 ] improves C3 plant growth rate and results in consequent increases in crop yield. CO 2 -enriched air may shift plant-pathogen interaction and plant disease resistance via several distinct mechanisms: alteration of plant primary metabolic status, up-or down-regulation of critical resistance genes, modification of stomatal opening or density and changes in pathogen virulence (reviewed by Mhamdi and Noctor 2016). In contrast to C3 plants, the growth and yield of C4 plants cannot be enhanced directly (Vaughan et al. 2014). It might be due to the presence of Kranz anatomy and unique CO 2 concentrating mechanism (CCM) (Gowik and Westhoff 2011). The primary metabolic pathways triggered by elevated [CO 2 ] can also act as a feedback to plant defence response (Mhamdi and Noctor 2016). Atmospheric [CO 2 ] increase and consequent global warming and drought also induce the production of secondary metabolites that are key to plant immune responses, i.e. salicylic acid (SA) in C3 plants (Casteel et al. 2012) and terpenoid phytoalexins in C4 plants (Vaughan et al. 2015). As reported, elevated [CO 2 ] makes maize, the most important C4 crop, more susceptive to fungi F. verticillioides. It is a consequence of the down-regulation of lipoxygenases (LOX) genes and attenuation of jasmonic acid (JA) and terpenoid phytoalexins production (Vaughan et al., 2014(Vaughan et al., , 2016. In recent years, high through-put transcriptome profiling has been widely used to provide first clues to the genetic changes in plants under elevated [CO 2 ] conditions. Microarray analysis was performed in Arabidopsis, poplar, soybean and wheat to investigate high [CO 2 ]-induced transcriptomic alterations. Prins et al. reported that 2× [CO 2 ] altered leaf transpiration rates, carbohydrate metabolism and protein carbonyl accumulation in maize in a leaf-rank specific manner by using a couple of transcriptome and proteome analysis (Prins et al. 2011). RNA sequencing (RNA-seq) technology is a robust and convenient platform that provides an unbiased view for the transcriptome of both model and non-model organisms (Wang et al. 2009). Here we employed RNA-seq to reveal the transcriptomic differences among maize leaf tissues under elevated [CO 2 ] condition for 0, 7 and 14 d. However, in nature, the adaptative response to the increasing atmospheric [CO 2 ] occurs during the whole life cycle in plants. We realized the limitation in short exposure time course in growth chambers since plants adaptation in natural environments may require even longer time. As reported, carbon sink rate in land plants is influenced by the atmospheric [CO 2 ] (Hendrey et al. 1993). Thus, we attempted to use a high CO 2 concentration (ten times) to accelerate the process. The discovery in this project provides a first glimpse of the genes that determine the high [CO 2 ]-induced adaptive response in maize and will also facilitate the research in molecular markers identification and genetics feeding in maize agriculture.

Plant materials and growth condition
The hybrid maize Zea mays seeds var. DH305 were purchased from Denghai Seeds Co. Ltd. (Yantai, China). The germination and growth conditions were 10 h light, 25°C/ 14 h dark, 22°C cycles in a growth chamber (chamber A). Medical-grade sterile air containing 4000 μmol CO 2 mol −1 was pre-aerated into another growth chamber (chamber B) constantly for at least 24 h before experiments. Eight-week-old maize seedlings was subsequently moved from chamber A to B. Chamber B was aerated with air containing 4000 ppm CO 2 for a further 14 d continuously. The CO 2 concentration in growth chamber B was monitored by an infrared CO 2 analysis reader (Huayun, China). Relative humidity of growth chamber was maintained at 30% ± 5% in the lab room with fixed air humidity. Maize leaf tissues were harvested after 0, 7 and 14 d in elevated [CO 2 ] chamber for analysis.

Total RNA isolation and cDNA library construction
Total RNA was isolated using RNeasy Kit for Plants (QIA-GEN, Germany). mRNA was subsequently enriched and fragmented. The fragmented mRNA was used as templates in cDNA first strand synthesis. The anti-sense cDNA strands were produced by the addition of dNTPs, RNase H and DNA polymerase I. The double-stranded cDNA was purified and end-repaired, adenylated and ligated to sequencing adaptors. Targeted fragments were reclaimed through agarose electrophoresis and amplified in PCR.

High throughput sequencing, data analysis and gene functional annotation
The cDNA libraries were sequenced using PE150 strategy on Illumina HiSeq 2000 platform (Illumina, USA), which generated a total of 16GB raw data. Raw data were further filtered to eliminate low-quality reads, adaptor polluted reads, and ambiguous Ns reads and to produce clean data. Once clean, the data were aligned to the maize Zea mays genome (NCBI BioProject Accession NO. PRJNA10769) and assembled using the software TopHat v2.0.12 ), which deployed BowTie2 results to align splice junctions (Langmead et al. 2009). The significantly DEGs were selected by using the fold change > 2 and averaged value > 1 of the normalized read counts of genes at the two-time points. GO terms and KEGG pathways enrichment analysis on the DEGs were performed using Metascape Arabidopsis thaliana database and ArgiGO Zea mays ssp. Database (Tian et al. 2017). The UniProt accession numbers of DEGs were used as input for the analysis using both Metascape and AgriGO.
Transcript levels were calculated as relative to the gene of interest/β-TUB4 ratio, and the value was subsequently presented as a relative percentage of control. Error bars indicate the ± S.D. value for triplicates. One-Way ANOVA was used to determine statistical differences (significant value at p < 0.05). Primers used in qRT-PCR analysis are summarized in Table 3.

Results
Maize seedlings in the high-CO 2 chamber were inspected for evidence of bacterial or fungi infection and in all tissues, we did not find evidence of pathogen invasion i.e. lesions. In comparison with seedlings in normal [CO 2 ] condition, maize seedlings conditioned under 10× [CO 2 ] stress did not exhibit notable difference in phenotype.

Degs highlights the effect of high [CO 2 ] on plant immunity
As reported, higher [CO 2 ] altered the expression of several JA biosynthesis-related genes i.e. lipoxygenases (LOXs), allene oxide synthases (AOSs), allene oxide cyclases (AOCs) and SA-mediated defense regulatory gene NONEXPRESSOR of pathogenesis-related gene1 (NPR1) in maize stem tissues and eventually increased maize susceptibility to fungal pathogens (Vaughan et al. 2014). Thus, we further analysed the expression level of seven DEGs: LOX2, LOX5, LOX12, AOS1, AOS2, AOC1 and NPR1 from three sets using qRT-PCR. As shown in Figure 1, qRT-PCR data are generally consistent with RNA-seq data, providing an independent confirmation of expression level alteration.

Elevated atmospheric [CO 2 ] initiates alterations in a couple of biological processes
The relative expression level of all DEGs at the three distinct time points is illustrated in Figure 2(a). The top GO term  obtained from Metascape Arabidopsis thaliana database in the three comparison sets is the same, namely 'protein phosphorylation' (Figure 2(b-d)). There are 90, 126, 132 up-regulated genes and 83, 131, 121 down-regulated genes in this process in the three comparison sets, respectively. Other enriched GO terms, such as 'protein ubiquitination', 'plant organ development' 'cellular response to endogenous stimulus' and 'cell death' were highlighted in all of the three comparison sets. 'Defense response to bacterium' was identified in both of the Zm_E1 (7 d) to Zm_C (0 d) and Zm_E2 (14 d) to Zm_C (0 d) comparison set. Additionally, the top GO term obtained from the GO terms enrichment analysis using AgriGO maize Zea mays ssp. database was 'oxidationreduction' (Figure 3), which was also reported in a recent study (Mhamdi and Noctor 2016). The top KEGG pathway enriched from the genes supporting the GO term 'protein phosphorylation' in the three comparison sets is 'plant-pathogen interaction' (Figure 4), which is consistent with the previously reported results that elevated CO 2 concentration has great impact on the defence system of maize (Vaughan et al. 2014). In the KEGG map of 'Plantpathogen interaction', we identified 33, 51 and 38 up-regulated genes and 13, 24, 26 down-regulated genes in each of the comparison set, respectively. The disease resistant protein, RPM1, RIM1-interacting protein 4 (RIN4) and RPS2 were up-regulated while LRR receptor-like serine/ threonine kinase proteins were down-regulated in this KEGG map, indicating that elevated CO 2 concentration stress may affect the plant hypersensitive reaction (HR).

Gene interaction networks participate in elevated [CO 2 ]-induced stress
Three densely connected networks showing several genes/proteins interactions were found in the DEGs using STRING protein interaction database (http://string-db.org) (Figure 5 (a-c)). As presented in Figure 5(d), the expression level of transcription factor c-Myc gene: MYC−2, −4 and three JA pathway repressor, Jasmonate-ZIM-domain protein genes: JAZ−1, −2 (TIFY10B) and −3 all increased at 7d, with a slight decrease at 14 d in comparison with control. On the other hand, JAZ −12 showed an apparent reduction at 7 d and a subsequent increase at 14 d. Figure 5(e) shows the relative expression level of all of the 17 WRKY genes in the DEGs sets, most of which have increased expression at 7 d in comparison with control and followed by the decrease in expression at 14 d. Other WRKY genes also changed with different change profiles in the three comparison sets. Furthermore, the up-regulation of WRKY−22 and −33 were identified in the KEGG map of 'plant-pathogen interaction'. As such, it probably suggests these two WRKY have different functions in plant resistance that remains poorly understood and are activated in response to elevated [CO 2 ] stress. Figure 5(f) shows the relative expression level of three TCP genes: TCP−5, −11 and −15, together with the transcription factor PYE.

Discussion and conclusion
Plants participate in the global carbon cycle through CO 2 uptake and carbon fixation and have emerged as predominant modulators of atmospheric [CO 2 ]. As reviewed above, alteration in atmospheric [CO 2 ] has influenced plants development, growth, immunity and crop harvests. Here we try to illustrate an in-depth view of maize leaf tissues transcriptome conditioned under elevated [CO 2 ] stress in order to analyse DEGs and to decipher hampered physiological perspectives. We selected a relatively high [CO 2 ] stress, ten times of ambient atmospheric CO 2 concentration, and a relatively longer exposure time (7 d and 14 d) in our experiments in comparison with previous reports. Some of the changes we found are  similar to the published data in C3 plants under moderate [CO 2 ], therefore validating our approach. Furthermore, novel discovery has been also provided from our results, which are now a basis for further experiments.
In contrast to previous reports, our RNA-seq data could record the changed profiles in the expression level of several thousand transcripts within 14 d. For the first time, the importance of two protein post-translational modification pathway, 'protein phosphorylation' and 'protein ubiquitination' is highlighted in GO enrichment results. The phosphorylation and ubiquitination could alter protein conformation and stability, which eventually determine enzyme activity and cellular homeostasis (Karve and Cheema 2011). Thus, maize proteome analysis using iTRAQ and phosphoproteome analysis using LC-MS/MS will provide more clues. Out of our KEGG pathway analysis, we verified the connection between elevated [CO 2 ] stress and plantpathogen interaction and MAP kinase signaling pathway.
C-Myc is a basic helix-loop-helix (bHLH) transcription factor and direct target of JA signaling repressor,  (Xu et al. 2006;Chen et al. 2010). The interaction among three WRKY genes may imply that defence to pathogens in maize is affected. The TEOSINTE BRANCHED1/ CYCLOIDEA/PCF (TCP) transcription factor family is responsible for triggering plant leaf maturation and plant immunity (Aguilar-Martinez and Sinha 2013). POPEYE (PYE) is another bHLH transcription factor that regulates several biotic and abiotic stresses, i.e. iron deficiency in plants (Long et al. 2010;Kim et al. 2014). The interaction within TCP and PYE may be engaged in the regulation of hormone CK and iron-deficiency related stimulus in plant organ development. Collectively, the three gene clusters are consistent with previous reports of the essential function of plant hormone production and plant defence response to pathogens in high [CO 2 ] condition.
Free-Air Carbon dioxide Enrichment (FACE) is a powerful method developed by Brookhaven National Laboratory for ecologists to measure plant growth under elevated atmospheric [CO 2 ] in a specified area outdoors (Hendrey et al. 1993). In comparison with greenhouse experiments, FACE experiments are close to natural conditions in the evaluation of increased [CO 2 ] effects and could take into account factors such as plant competition, that are missing from lab experiments. However, FACE experiments require longer time (one or two years) and larger field area (Dier et al. 2018) than lab experiments. Our experiments in growth chamber is different from FACE experiments but have allowed us to pinpoint areas than can be now investigated using a FACE design with more realistic conditions.
In conclusion, the utilization of 10× [CO 2 ] stress and relatively prolonged exposure time courses in this project enable us to identify 1390 DEGs and three transcription factors clusters. The supporting data provided in supplementary materials with this article can also serve as a genetic database for scientists to screen genes of interest.