Transcriptional Analysis of Masson Pine (Pinus massoniana) under High CO2 Stress

To explore the molecular mechanism of the response of Masson pine (Pinus massoniana), the main coniferous tree in southern China, to high CO2 stress, transcriptome sequencing was carried out to analyze the genome-wide responses of annual seedlings under different durations (0 h, 6 h, 12 h and 24 h) of high CO2 stress. The results showed that a total of 3080/1908, 3110/2115 and 2684/1483 genes were up-/down-regulated after 6 h, 12 h and 24 h of treatment, respectively, compared with control check group (CK, 0 h). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that most of these differentially expressed genes (DEGs) were enriched in energy metabolism, carbohydrate synthesis, cell wall precursor synthesis and hormone regulation pathways. For energy metabolism, the expression of most genes involved in photosynthesis (including the light reaction and Calvin cycle) was generally inhibited, while the expression of genes related glycolysis, the tricarboxylic acid (TCA) cycle and PPP pathway was up-regulated. In addition, the increase in the CO2 concentration induced the up-regulation of gene expression in the sucrose synthesis pathway. Among all starch synthesis genes, GBSS (granule-bound starch synthase) had the highest expression level. On the other hand, during the synthesis of hemicellulose and pectin (cell wall precursor substances), the expression levels of GMD (GDP-mannose 4,6-dehydratase), MGP (Mannose-1-phosphate guanylyl transferase) and RHM (Rhamnose biosynthetic enzyme) were the highest, suggesting that the synthesis of the raw materials hemicellulose and pectin in Masson pine under stress were mainly supplied by GDP-Man, GDP-Fuc and UDP-Rha. Finally, stress inhibited gene expression in the ABA (Abscisic Acid) synthesis pathway and induced gene expression in the GA (Gibberellin), SA (Salicylic acid), BR(Brassinolide) and MeJA (Methyl Jasmonate) pathways. Stomatal switches were regulated by hormonal interactions. This experiment elaborated on the response and molecular mechanism of Masson pine to CO2 stress and aided in screening carbon sequestration genes for the corresponding molecular research of Masson pine in the future.


Introduction
As humans enter an industrial society, environmental issues have become increasingly prominent. Global warming has been widely monitored by governments and the public worldwide and has risen to become one of the most important political, diplomatic, and economic issues [1]. According to the Global Carbon Project (GCP), global CO 2 emissions from fossil fuels and industry reached 36.8 Gt in 2017, an increase of approximately 65% from the baseline year (1990) of the Kyoto Protocol [2].

Plant Material and Experimental Conditions
One-year-old Masson pine seedlings, obtained from the seed orchard of Baisha state-owned forest farm, Shanghang, Fujian Province, China (25 • 15 N, 116 • 62 E), were used in this study. Individuals of the same clones with similar heights, uniform growth and strong growth potential were selected as the test materials and subsequently moved into a growth chamber to recover for 15d. The growth conditions were 10 h light/14 h dark cycles at 25 • C in the chamber. Air containing about two times of the chamber CO 2 concentration before experiment was aerated into the growth chamber constantly for at least 24 h. The CO 2 concentration in the chamber was monitored by an infrared CO 2 analysis reader (SenseAir, Delsbo, Sweden). The seedlings were sampled after 6 h, 12 h and 24 h of treatment with the high CO 2 concentration, and the needles were selected for downstream experiments.

Total RNA Isolation, Complementary DNA Library Preparation and Sequencing
Total RNA from four treatments (0 h, 6 h, 12 h and 24 h) seedlings (Among them, 0h treatment was considered as the control check group (CK group)) with three biological replicates for each treatment was extracted using the Plant RNA Isolation Kit (Tiangen, Beijing, China). Sequencing library were constructed using RNA Library Prep Kit for Illumina (NEB, Boston, MA, USA). Then the libraries were sequenced using a Hiseq 4000 (Illumina, San Diego, CA, USA), and generated 150 bp paired-end reads. To get clean reads, sequences with length less than 30 bp, reads with N ratio over 10% and without inserted fragments due to reasons such as connector self-connection and adapter sequences were removed using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) [18]

Transcriptome Assembly and Functional Annotation
De novo assembly of all the clean reads was conducted using Trinity version 2.5.0 (https://github. com/trinityrnaseq/trinityrnaseq/wiki) based on all parameters set as their defaults [19]. Transcripts corresponding to paralogous genes were sorted out to finally obtain the assembly sequences. TransRate software version 1.0.3 (http://hibberdlab.com/transrate/) was used to filter and optimize the initial assembly sequences obtained from Trinity [20], and the assembled sequence was evaluated using BUSCO (Benchmarking Universal Single-Copy Orthologs) version 3.0.2 (http://busco.ezlab.org) [21]. Both of the two softwares ran with their default parameters.
To functionally annotate genes, the sequences were BLASTed in Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/) public databases. KEGG enrichment analysis for the DEGs was carried out by KOBAS version 3.0. In this software, False-positive were assessed by BH (FDR correction with Benjamini/Hochberg) method with a cut-off E-value of 10 −6 . When the p-adjust (the adjusted p-value) of a KEGG pathway was less than 0.05, we considered this KEGG pathway was significantly enriched.

Differential Expression Analysis
The gene expression level for each sample was determined according to the transcripts per million reads (TPM) using RSEM version 1.2.31 (Univ Wisconsin, Madison, WI, USA) with all parameters set as its defaults [22], in which the read counts were normalized using the edgeR package with the Trimmed Mean of M-values method, and then the length of the gene (L) and the normalized counts (N) were used to calculate the TPM (TPM = 10 6 × (N i /L i )/ (N j /L j )). The DEGseq R package version 1.10.1 was used to analyze differential expression of different samples [23]. The significant differential expression threshold was set as q-value < 0.005 and |log 2 (foldchange)| ≥ 1 [24]. The heatmap and differential expression of genes among samples were analyzed using "pheatmap" R package version 1.0.12 (Massachusetts General Hospital, Boston, MA, USA).

Quantitative Real-Time PCR Validation
To validate the RNA-sequencing (RNA-seq) results, the expression levels of nine genes were measured using qRT-PCR. Reaction mixtures consisted of 10 µL of 2× TOP Green qRT-PCR SuperMix (TOYOBO Biotech, Shanghai, China), 0.4 µL of forward primer and reverse primer, 2 µL of complementary DNA (cDNA), 0.4 µL of 50× Passive Reference Dye and 6.8 µL of ddH 2 O. The PCR program was set up in six stages: (1) 94 • C for 30 s (preincubation), (2) 94 • C for 5 s, (3) 55 • C for 15 s, (4) 72 • C for 10 s, repeated 40 times (amplification), (5) 95 • C for 0.5 s, and (6) 60 • C for 1 min (melt). The PCR quality was estimated based on melting curves. TUA (Alpha-tubulin) was used as the internal control [25]. The gene-specific primers employed are shown in Table A1 from Appendix A. Three independent biological replicates and three technical replicates for each biological replicate were run. Quantification was achieved using comparative cycle threshold (Ct) values, and gene expression levels were calculated using the 2 −∆∆Ct method [25]. The significance was determined by t-test using SPSS statistical software (IBM, New York, NY, USA) (p < 0.05).

Transcriptome Sequencing and De Novo Assembly
The cDNA libraries of the four treatments (0 h (CK), 6 h, 12, and 24 h) were sequenced and generated a total number of 49,314,299, 47,459,322, 45, The Trinity software generated 140,863 transcripts with an average length of 891 bp and an N50 of 1463 bp. In total, 92,424 unigenes were obtained in the range of 201~15,491 bp. Of these, 48,592 (52.57%) were less than 500 bp, 22,267 (24.09%) were 501~1000 bp, 12,887 (13.94%) were 1001~2000 bp and the remaining 8678 (9.39%) were >2000 bp ( Table 2). TransRate and software BUSCO evaluated the assembly results, and transcript score was 0.20045 and 77.7%, respectively, and unigenes was 0.30498 and 74.2%, respectively ( Table 2). According to previous reports [26,27], in Masson pine, unigene obtained by transcriptome sequencing under other stress treatments ranged from 70,896 to 101,806. Our results were similar to them. Combined with the N50 and Q30, we believe that the sequencing results are relatively reliable and could be further analyzed.

Gene Expression and KEGG Enrichment Analysis under CO 2 Stress
A total of 7088 genes were differentially expressed between the samples from the three CO 2 stress treatments and the control samples. Of these, 4988, 5225 and 4167 genes were differentially expressed between 6 h and CK, 12 h and CK, and 24 h and CK treatments, respectively ( Figure 1A). Among the differentially expressed genes (DEGs), 3080/1908, 3110/2115 and 2684/1483 genes were up-/down-regulated at 6 h, 12 h, and 24 h, respectively, compared with CK ( Figure 1B). Gene enrichment analysis of the DEGs based on KEGG analysis revealed that these genes were mainly involved in several pathways at different time points, including photosynthesis, carbon fixation (the Calvin cycle), glycolysis, the tricarboxylic acid (TCA) cycle, starch and sucrose metabolism, fructose and mannose metabolism, galactose metabolism and plant hormone signal transduction, etc. ( Figure 1C,D,E). In the three KEGG enrichment Figures that compare the different treatment time points with CK, all of the above pathways ranked within the top 20 pathways (except the "TCA cycle" in 6 h/CK, Figure 1C), indicating that energy metabolism, carbohydrate synthesis, cell wall synthesis and hormone regulation may be the main metabolic activities in Masson pine under high CO 2 stress.

Gene Expression and KEGG Enrichment Analysis under CO2 Stress
A total of 7088 genes were differentially expressed between the samples from the three CO2 stress treatments and the control samples. Of these, 4988, 5225 and 4167 genes were differentially expressed between 6 h and CK, 12 h and CK, and 24 h and CK treatments, respectively ( Figure 1A). Among the differentially expressed genes (DEGs), 3080/1908, 3110/2115 and 2684/1483 genes were up-/down-regulated at 6 h, 12 h, and 24 h, respectively, compared with CK ( Figure 1B). Gene enrichment analysis of the DEGs based on KEGG analysis revealed that these genes were mainly involved in several pathways at different time points, including photosynthesis, carbon fixation (the Calvin cycle), glycolysis, the tricarboxylic acid (TCA) cycle, starch and sucrose metabolism, fructose and mannose metabolism, galactose metabolism and plant hormone signal transduction, etc. ( Figure  1C,D and E). In the three KEGG enrichment Figures that compare the different treatment time points with CK, all of the above pathways ranked within the top 20 pathways (except the "TCA cycle" in 6 h/CK, Figure 1C), indicating that energy metabolism, carbohydrate synthesis, cell wall synthesis and hormone regulation may be the main metabolic activities in Masson pine under high CO2 stress.  Figure 2 shows the energy metabolic pathways and the general pattern of the relative changes in the related gene expression patterns in Masson pine under high CO 2 concentration conditions. Notably, in the transcriptome data, the pathways involved in energy metabolism were polarized. After Masson pine had been exposed to CO 2 stress, in the photosynthesis pathways, including the photoreactions and Calvin cycle (Figure 2A, black arrow), the relative expression of each gene in the metabolic pathway except for GAPD (glyceraldehyde-3-phosphate dehydrogenase) showed a downward trend ( Figure 2B). GAPD catalyzes 3-phosphoglyceraldehyde (PGAL) to 1,3-diphosphoglycerate (1,3-DPG) in the Calvin cycle and 1,3-DPG to dihydroxy-acetone phosphate (DHAP) in the Embden-Meyerhof-Parnas (EMP) pathway ( Figure 2A, orange arrow). Since GAPD plays a role in these two metabolic pathways, the change in its expression patterns may be the result of superposition. Moreover, the gene expression patterns involved in both photosynthesis and other metabolic pathways (such as RPI (Ribulose Phosphate Isomerase) or RPE (Ribulose Phosphate Epimerase)) decreased more slowly than those involved only in photosynthesis (such as RCA (Rubisco Activase) or FBP (Fructose-1,6-diphosphate)). On the other hand, the relative expression levels of genes in other pathways involved in energy anabolism, including the EMP pathway, the pentose phosphate (PPP) pathway ( Figure 2A, blue arrow) and TCA (Figure 2A, green arrow), generally increased ( Figure 2B). Among them, the expression rate of OGDC (α-ketoglutarate dehydrogenase) increased the fastest and was 4.8, 6.09 and 6.18 times that of CK at 6 h, 12 h and 24 h, respectively. In general, the omics data revealed that energy metabolism was strongly enhanced for contributing to elevated CO 2 tolerance in Masson pine, except for the light reaction and Calvin cycle.

Energy Metabolism under Elevated CO 2 Stress
photoreactions and Calvin cycle (Figure 2A, black arrow), the relative expression of each gene in the metabolic pathway except for GAPD (glyceraldehyde-3-phosphate dehydrogenase) showed a downward trend ( Figure 2B). GAPD catalyzes 3-phosphoglyceraldehyde (PGAL) to 1,3-diphosphoglycerate (1,3-DPG) in the Calvin cycle and 1,3-DPG to dihydroxy-acetone phosphate (DHAP) in the Embden-Meyerhof-Parnas (EMP) pathway ( Figure 2A, orange arrow). Since GAPD plays a role in these two metabolic pathways, the change in its expression patterns may be the result of superposition. Moreover, the gene expression patterns involved in both photosynthesis and other metabolic pathways (such as RPI (Ribulose Phosphate Isomerase) or RPE (Ribulose Phosphate Epimerase)) decreased more slowly than those involved only in photosynthesis (such as RCA (Rubisco Activase) or FBP (Fructose-1,6-diphosphate)). On the other hand, the relative expression levels of genes in other pathways involved in energy anabolism, including the EMP pathway, the pentose phosphate (PPP) pathway (Figure 2A, blue arrow) and TCA (Figure 2A, green arrow), generally increased ( Figure 2B). Among them, the expression rate of OGDC (α-ketoglutarate dehydrogenase) increased the fastest and was 4.8, 6.09 and 6.18 times that of CK at 6 h, 12 h and 24 h, respectively. In general, the omics data revealed that energy metabolism was strongly enhanced for contributing to elevated CO2 tolerance in Masson pine, except for the light reaction and Calvin cycle.   Table A2 from Appendix A.

Biosynthesis of Sucrose, Starch and Cell Wall Components under Elevated CO 2 Stress
The expression levels of the genes in the sucrose and starch synthesis pathways generally showed a trend of up-or slightly down-regulation under CO 2 stress conditions compared with CK, except for AGP (Adenosine Diphosphoglucose Pyrophosphorylase) and SBE (Starch Branching Enzyme) ( Figure 3B). AGP catalyzes adenosine diphosphate glucose (ADP-Glc) from glucose-1-phosphate (G-1-P), and its expression continued to decline with increasing stress time. On the other hand, among the genes in starch metabolic pathway, the expression of GBSS (Granule-bound Starch Synthase) was higher than SSS (Soluble Starch Synthase) and SBE.
After a series of reactions, triose phosphate was transformed into UDP-Glc in the cytoplasm ( Figure 3A), and then it was catalyzed to sucrose by sucrose-phosphate synthase (SPS). After being subjected to stress treatment, the expression levels of SPS at 6 h, 12 h and 24 h were no significant difference ( Figure 3B). Meanwhile, sucrose can be transformed to UDP-Glc and hexose by sucrose synthase (SUS) and cell wall invertase (cwINV), and hexose can then be catalyzed by cytosolic invertase (INV) to form fructose-6-phosphate (F-6-P) and glucose-6-phosphate (G-6-P), which are the synthetic precursors of UDP-Glc ( Figure 3A). In the above series of reactions, the expression of SUS, cwINV and INV were up-regulated in different degrees under CO 2 stress compared with CK. In addition, UDP-Glc is the precursor for cellulose synthesis, which is catalyzed by cellulose synthase complex (CSC), including cellulose synthase subunit (CesA), cellulose synthase (Csl) and its isoenzymes ( Figure 3A). The experimental results showed that the expression levels of CesA, Csl A, Csl C and Csl D were up-regulated with increasing stress time. However, Csl B and Csl E decreased at the same time ( Figure 3B). In general, the omics data revealed that gene expression patterns in sugar and cell wall component metabolic pathways were enhanced under elevated CO 2 stress in Masson pine.
After a series of reactions, triose phosphate was transformed into UDP-Glc in the cytoplasm ( Figure 3A), and then it was catalyzed to sucrose by sucrose-phosphate synthase (SPS). After being subjected to stress treatment, the expression levels of SPS at 6 h, 12 h and 24 h were no significant difference ( Figure 3B). Meanwhile, sucrose can be transformed to UDP-Glc and hexose by sucrose synthase (SUS) and cell wall invertase (cwINV), and hexose can then be catalyzed by cytosolic invertase (INV) to form fructose-6-phosphate (F-6-P) and glucose-6-phosphate (G-6-P), which are the synthetic precursors of UDP-Glc ( Figure 3A). In the above series of reactions, the expression of SUS, cwINV and INV were up-regulated in different degrees under CO2 stress compared with CK. In addition, UDP-Glc is the precursor for cellulose synthesis, which is catalyzed by cellulose synthase complex (CSC), including cellulose synthase subunit (CesA), cellulose synthase (Csl) and its isoenzymes ( Figure 3A). The experimental results showed that the expression levels of CesA, Csl A, Csl C and Csl D were up-regulated with increasing stress time. However, Csl B and Csl E decreased at the same time ( Figure 3B). In general, the omics data revealed that gene expression patterns in sugar and cell wall component metabolic pathways were enhanced under elevated CO2 stress in Masson pine.   [28] and Jana [29]. The blue, purple and green rectangles represent the cytosol, chloroplast and Golgi, respectively. The different enzymes are shown in red font. (B) Expression changes in the genes involved in metabolic pathways in response to stress. White indicates no change, red up-regulation, and blue down-regulation in each treatment, as shown in the color bar for a log 2 fold change scale. The abbreviations in the figure are shown in Table A2 from Appendix A.

Hormone Regulation under Elevated CO 2 Stress
Based on previous studies [30][31][32][33], combined with KEGG analysis results, the expression patterns of 5 plant hormones (ABA, GA, SA, BR and JA) and their corresponding synthesis genes under CO 2 stress were analyzed in this study. In the ABA synthesis pathway, except for the expression of ZEP (zeaxanthin epoxidase) that was slightly up-regulated (no significant difference from CK), all other genes were down-regulated. Among them, NCED (9-cis-epoxycarotenoid dioxygenase) decreased most obviously ( Figure 4B). and the log2 fold changes in expression compared to CK decreased 4.73, 6.10 and 5.38 times at 6 h, 12 h and 24 h, respectively.
On the other hand, the genes expression pattern of the other four hormones showed an increasing trend under CO 2 stress. The key enzyme in the GA metabolic pathway is geranyl geranyl pyrophosphate (GGPS), which catalyzes the synthesis of geranylgeranyl pyrophosphate (GGPP) from isopentenyl pyrophosphate (IPP), geranyl pyrophosphate (GPP) and farnesyl pyrophosphate (FPP). Under CO 2 stress, GGPS expression increased continuously with time ( Figure 4B).
As shown in Figure 4A, SA can be synthesized via two routes, by isochorismic acid (ICA) and by cinnamic acid (Ca), benzoic acid (Ba), etc. The specific genes involved in the first pathway have not yet been thoroughly studied [34], while the second pathway has been well researched. Under CO 2 stress, PAL expression continued to increase, reaching a peak at 12 h and then stabilizing (there was no significant difference between 24 h and 12 h). Due to the expression pattern of PAL, we speculated that the pathway used for SA synthesis in Masson pine under CO 2 stress was mainly the second pathway.
The BR synthesis pathways include early ( Figure 4A, yellow block) and advanced (blue block) C-6 synthetic pathways [35]. Intermediate metabolites in the advanced pathway could be converted to corresponding metabolites in the early pathway. However, at the upstream pathway, the conversion efficiency was not very high because the expression levels of ROT (C-23 hydroxylase) was down-regulated with increasing treatment time. In contrast, downstream of the pathway, due to the increase in CYP expression ( Figure 4B), BR might be synthesized through the early pathway. Overall, both metabolic pathways showed increased efficiency under CO 2 stress ( Figure 4B).
All genes in the JA (MeJA) synthetic pathway were up-regulated under CO 2 stress. Among them, AOC (allene oxide synthase) was the most significantly up-regulated. At 6 h, the relative expression level of AOC was 2.84 times that of CK and then increased, and may continue to increase with increased treatment time. Within 24 h, the expression of AOC was far higher than that of the other genes in the pathway ( Figure 4B). It could be speculated that AOC might be a key gene involved in the induction of JA production in Masson pine under high CO 2 concentrations.   Table A2 from appendix A.

Validation by Quantitative Real-Time PCR
To verify the reliability of the transcriptome data, nine genes showing significant up-or down-regulation in the stressed seedlings were randomly chosen for qRT-PCR analysis ( Figure 5). Among them, PAL showed constitutively up-regulated expression, and CYP increased at first and reached maximum expression at 6 h, after which there was no significant change. Compared with The hormone biosynthesis pathways and stomatal regulation mechanism according to Zhao [35]. The mazarine, green, brown, purple and wathet frame represents ABA, GA, SA, JA and BR synthetic pathways, respectively. The blue and yellow block represents advanced and early C-6 oxidation pathway in BR biosynthesis, respectively. The different enzymes are shown in red font. Sharp and T-shaped arrows indicate positive and negative regulation, respectively. (B) Expression changes in the genes involved in metabolic pathways in response to stress. White indicates no change, red up-regulation, and blue down-regulation in each treatment, as shown in the color bar for a log 2 fold change scale. The abbreviations in the figure are shown in Table A2 from Appendix A.

Validation by Quantitative Real-Time PCR
To verify the reliability of the transcriptome data, nine genes showing significant up-or down-regulation in the stressed seedlings were randomly chosen for qRT-PCR analysis ( Figure 5). Among them, PAL showed constitutively up-regulated expression, and CYP increased at first and reached maximum expression at 6 h, after which there was no significant change. Compared with CK, 6 h, 12 h and 24 h CO 2 treatments had up-regulated levels of RHM and OPR, and the relative expression levels of these genes were higher at 6 h than at 12 h and 24 h. On the other hand, enolase had the same trend as RHM and OPR, except that the maximum value appeared at 12 h. CPS and SBP showed constitutively down-regulated expression, and TPI (triose-phosphate isomerase) showed a trend of first decline and then increase. Among all genes, the qRT-PCR results of 6 genes (CYP, CPS, OPR, PAL, SBP and TPI) were very close to the RNA-seq (via the TPM algorithm) results. The expression trend of RHM and enolase was similar to that of RNA-seq, but the fold change in qRT-PCR expression was lower than that in the RNA-seq data. The qRT-PCR results of CesA were quite different compared with the RNA-seq results. In general, the RNA-seq data were similar to the gene expression trend shown by qRT-PCR analysis, indicating that the results of RNA-seq analysis were effective.

Effects of Elevated CO2 on Masson Pine Energy Metabolism
The results of this study showed that in the major energy metabolism pathways of Masson pine, the expression levels of most genes involved in photosynthesis, including light and dark reactions, showed a significant downward trend under CO2 stress. We speculate the main cause of this phenomenon should be the decrease in the expression of Rubisco, which is a key enzyme in photosynthetic carbon assimilation, directly affecting the photosynthetic efficiency of Masson pine [36]. Combined with the up-regulated of key genes in the carbohydrate metabolism during the same period ( Figure 3B), it leads to the activation of its signaling pathway, which inhibits the expression Figure 5. The expression changes in the 9 randomly selected genes were determined using quantitative real-time PCR (qRT-PCR) results and sequencing data. The x-axis represents different processing times, and the y-axis represents changes in gene expression under CO 2 stress. The data show the fold change in the expression of each gene under high CO 2 relative to control conditions. Error bars represent standard deviations. Red indicates the RNA-sequencing results under the TPM (transcripts per million reads) algorithm, and blue indicates the qRT-PCR results.

Effects of Elevated CO 2 on Masson Pine Energy Metabolism
The results of this study showed that in the major energy metabolism pathways of Masson pine, the expression levels of most genes involved in photosynthesis, including light and dark reactions, showed a significant downward trend under CO 2 stress. We speculate the main cause of this phenomenon should be the decrease in the expression of Rubisco, which is a key enzyme in photosynthetic carbon assimilation, directly affecting the photosynthetic efficiency of Masson pine [36]. Combined with the up-regulated of key genes in the carbohydrate metabolism during the same period ( Figure 3B), it leads to the activation of its signaling pathway, which inhibits the expression of RbcS (Rubisco small subunit) and other photosynthesis-related protein-coding genes through hexokinase transmission, thereby affecting the photosynthetic efficiency [37].
Significantly, previous studies [38,39] have shown the expression of Rubisco will show an up-regulated trend under the high CO 2 treatment condition. But the results of this experiment are the opposite. This may be because the expression of carbon sequestration gene is greatly affected by the leaf age [40]. The photosynthetic acclimation of coniferous leaves in one-year needles under high CO 2 treatment was more obvious than that in mature coniferous leaves [41]. According to the experimental results, the expression of each gene in photosynthesis showed a downward trend when the sample was analyzed at the shortest detection time (6 h). Therefore, we speculate the photosynthetic acclimation should be completed within 6 h.
However, in other energy metabolism pathways of Masson pine, the expression levels of various genes under high CO 2 stress generally showed an upward trend. Fukayama et al. [42] obtained similar results in an experiment on the effect of CO 2 on rice, but the difference in gene expression was inferior to our study, possibly because rice belongs to herbage, and the adaptability of rice to CO 2 stress is higher than that of Masson pine. Previous studies found that a high concentration of CO 2 promoted the production of a large amount of sucrose in leaves, and the increase in sucrose concentration induces the expansion of glycolytic flux in plants and simultaneously increases the distribution of carbon in the organic acid and TCA cycles, thereby regulating the increase in related gene expression [43].

Effects of Elevated CO 2 on Masson Pine Carbohydrate and Cell Wall Component Synthesis
In photosynthetic cells, fixed carbon eventually transforms into sucrose and starch. These compounds are the main form of carbohydrate transport and storage in advanced plants. SPS, SUS and INV are considered to be the key genes in the sucrose synthesis pathway [44]. The results of our experiment suggest that there were no significant differences in SPS and SUS levels between the CK and different treatments, while INV showed an increased difference. This result suggested that, sucrose produced by photosynthesis was more likely to be decomposed into glucose-6-phosphate (G-6-P) and fructose-6-phosphate (F-6-P), and then participated in subsequent reactions is the main metabolic direction of sucrose under high CO 2 stress. Van et al. [45] found that INV was significantly up-regulated and SUS was relatively inhibited in mature leaves of in vitro tomato plantlets with high CO 2 stress under the condition of providing exogenous sucrose (3%). In addition, other studies also showed that the expression of INV and the contents of glucose and fructose in leaves increased under cold, NaCl, PEG and other stress conditions, which showed similar phenomena to CO 2 stress conditions [46,47].
In the starch synthesis pathway, three different starch-formed synthetases (GBSS, SSS and SBE) are the key regulatory enzymes. In this study, the expression level of GBSS were significantly higher than that of SBE and SSS under high CO 2 stress. GBSS generates a branchless linear glucan starch chain (amylose) by specific binding to starch granules, which is necessary for the production of amylose [48]. Due to the high expression of GBSS, we suspected that, amylose was more inclined to synthesize under high CO 2 stress in Masson pine. The similar results were confirmed in Cui [49] study. He found that the CO 2 stress could significantly improve the activity of GBSS in winter wheat, resulting in a large amount of amylose production, and the expression level increased first and then decreased, which was similar to this experiment.

Effects of Elevated CO 2 on Masson Pine Cell Wall Component Synthesis
The plant cell wall is mainly composed of polysaccharides, which are the largest storage place for photosynthetic carbon fixation [50]. The synthesis of cell wall components is a highly complex process involving multiple enzymes and metabolic intermediates [29]. The major components of the cell wall, such as cellulose, hemicellulose, and pectin, require a variety of ribose compounds, most of which rely on UDP-Glc derivatization ( Figure 3A) [29,51]. For the genes regulating UDP-Glc to produce hemicellulose and pectin precursor derivatives, the expression levels of GMD and RHM were higher than those of others. Due to the increased expression of RHM and GMD under CO 2 stress, we speculate that their corresponding products (UDP-Rha and GDP-Fuc) were the main components of cell wall precursors at this stage. However, there is still no unified conclusion on its specific process of influence. Sufficient long-term research and evidence to fully and thoroughly understand the mechanism are required.

Effects of Elevated CO 2 on Masson Pine Hormones and Stomatal Regulation
In this study, the results showed that except ZEP, the expression of each gene in the ABA synthesis pathway tended to decrease with prolonged CO 2 treatment time, especially NCED. The catalysis of 9-cis-neoxanthin to Xanthoxin by NCED is a key step in the ABA metabolic pathway in plastids [52].Therefore, a decrease in the expression of NCED seriously affects the expression of other genes in the ABA metabolic pathway, which in turn affects ABA accumulation in the cell [53]. Previous studies have shown that ABA is greatly affected by BR when regulating stomatal opening and guard cell physiological states [54]. Elevated CO 2 led to a sharp increase in the expression of key genes in the BR metabolic pathway, such as CYP, which would inhibite ABA from binding to ABI (a kind of serine/threonine protein phosphatase), thereby weakening the effect and signal transmission of ABA [55].
For other hormones, the key genes in their synthetic pathway increased to varying degrees at the same stage. Among them, the expression of GGPS, the key gene of GA synthesis, continued to up-regulate with prolonged stress. Studies have shown that GA can divide plant hypocotyl epidermal cells, promote stomatal formation, regulate stomatal density [56], so as to maintain stomatal conductance and the transpiration rate under the condition of elevated CO 2 concentration. On the other hand, according to our result, the key gene in the JA and SA synthetic pathway, including AOC, LOX and PAL, remained at a high expression level with prolonged stress. It indicated that the reaction proceeded in the direction of JA and SA synthesis, which would cause stomatal closure [57][58][59].
In general, although previous studies have shown that ABA is the most important hormone to control stomatal switching in plants, in Masson pine under CO 2 stress, the expression level of each gene in the ABA synthesis pathway was basically inhibited. However, this did not affect stomatal closure in Masson pine because the genes regulated other hormones that promote stomatal formation or induce its closure were general up-regulated in the same environment.

Conclusions
The effect of rising CO 2 concentrations on plants is known, however, little research has been done on Masson pine. In this study, we tried to explore the molecular response of Masson pine under high CO 2 stress. The results showed that, the genes expression would generally decrease in photosynthesis pathway (light reaction and Calvin cycle), and generally increase in other energy metabolic pathways, including TCA, EMP and PPP. At the same time, Increased CO 2 concentration could also promote the gene expression in cell wall precursor synthesis pathway. In addition, CO 2 stress inhibited the genes expression in the ABA synthesis pathway, but increased in other hormones synthesis pathway (including BR, GA, SA and MeJA), which might regulate stomatal density and stomatal closure. As the first report on the high-throughput sequencing of CO 2 tolerant Masson pine, this study should provide novel insights into CO 2 tolerance genes and be a valuable molecular basis for study in Masson pine.

Acknowledgments:
We would like to thank Major biotechnology corporation (Shanghai, China) for assistance with sequencing services.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Appendix A Table A1. Primers used in this study.