Comprehensive Transcriptome Analysis of Rare Carpinus putoensis Plants under NO2 stress

We evaluated a transcriptome using high-throughput Illumina HiSeq sequencing and related it to the morphology, leaf anatomy, and physiological parameters of Carpinus putoensis putoensis under NO2 stress. The molecular mechanism of the C. putoensis NO2 stress response was evaluated using sequencing data. NO2 stress adversely affected the morphology, leaf anatomy, and total peroxidase (POD) activity. Through RNA-seq analysis, we used NCBI to compare the transcripts with nine databases and obtained their functional annotations. We annotated up to 2255 million clean Illumina paired-end RNA-seq reads, and 250,200 unigene sequences were assembled based on the resulting transcriptome data. More than 89% of the C. putoensis transcripts were functionally annotated. Under NO2 stress, 1119 genes were upregulated and 1240 were downregulated. According to the KEGG pathway and GO analyses, photosynthesis, chloroplasts, plastids, and the stimulus response are related to NO2 stress. Additionally, NO2 stress changed the expression of POD families, and the HPL2, HPL1, and POD genes exhibited high expression. The transcriptome analysis of C. putoensis leaves under NO2 stress supplies a reference for studying the molecular mechanism of C. putoensis resistance to NO2 stress. The given transcriptome data represent a valuable resource for studies on plant genes, which will contribute towards genome annotations during future genome projects.


Introduction
Nitrogen dioxide (NO 2 ) is a product of nitric acid, which is used in industrial manufacturing; millions of tons of NO 2 are produced each year [1]. At high temperatures, NO 2 is a maroon gas with a typically harsh odor, and it is a key contributor to air pollution [2]. NO 2 is also an important component of acid rain [3]. Its corrosivity and highly oxidative nature make it harmful to plant biochemical and physiological processes after entering plants through the stomata [4]. In wild environments, the ambient NO 2 level that wild plants might encounter is 180 ppb. Currently, there are two theories regarding the effect of NO 2 on plants. The first is that NO 2 can form plant organic nitrogen compounds by being metabolized and amalgamated in the nitrate assimilation pathway [5]. Approximately 33% of NO 2 -derived N (NO 2 -N) taken up by plants was modified into a previously unknown Kjeldahl-unrecoverable organic nitrogen (unidentified nitrogen) [6], which can be incorporated into the α-amino group of soluble free amino acids [7,8], thereby not causing harm to the leaves [9,10]. Mayer et al. [11] investigated the changes in the physiological functions of NO 2 at a 10 µL L −1 concentration in Arabidopsis (Arabidopsis thaliana) cells and found that 1 h NO 2 fumigation induced pathogen resistance in the plant [11]. The second theory is that the majority of plants have a low absorption capacity for NO 2 [12]. Although most studies have investigated the amino acid response after NO 2 stress, there are no known reports on gene expression responses to NO 2 stress.

-N incorporation into the
Carpinus putoensis is a species in the Betulaceae family measuring approximately 15 m (49 feet) tall. It survives as a single tree on Putuo Island on the Zhoushan archipelago in China. It is monoecious but still able to reproduce sexually in nature [13]. The Zhejiang Forestry Science Research Institute has researched the cultivation and breeding of C. putoensis [14]; although the seed characteristics of C. putoensis were investigated previously, those studies stressed the characterization of the complete chloroplast genome and nuclear ribosomal sequence data [15]. It is vital to study C. putoensis resistance to NO 2 exposure to conserve this endangered species and improve its tolerance for future applications as a novel road greening and ornamental plant. Therefore, in a previous study, we evaluated the photosynthesis and Chl fluorescence responses of C. putoensis leaves to different NO 2 (6 µL/L) exposure times, both in terms of leaf gas exchange and the functionality of photosynthetic measurements [16]. Additionally, the chlorophyll content, the behavior of the stomata, and the ultrastructure of chloroplasts were analyzed together to find potential relationships between the photosynthesis in the leaves and cell transformation under NO 2 stress. However, a relationship between the leaf anatomy and transcription in C. putoensis under NO 2 stress has not previously been reported.
Therefore, in the current research, we evaluated the leaf anatomy and transcriptome gene expression of C. putoensis leaves under NO 2 stress. The purpose of this study is to provide a theoretical reference on the effects of traffic pollution on green plants.

Plant Material and Growth Conditions
One-year-old C. putoensis seedlings were grown in pots measuring 30 cm (open top) × 15 cm (height) × 20 cm (flat bottom) that were filled with well-mixed vermiculite, peat, and garden soil (1:1:1, v/v/v). In accordance with the water evaporation rate of the soil described by Allen et al. [17], they were watered with tap water every three days, and 1 L of full-strength Hoagland nutrient solution at was used biweekly for seedling cultivation. Before NO 2 treatment, the plants were allowed to grow naturally for 2 months [16].

NO 2 Fumigation
Fumigation was performed according to the method described in the literature [11]. open-top NO 2 fumigation glass chambers measuring 50 × 50 × 50 cm were built. The plants were fumigated with NO 2 at 6 µL/L that was supplied by cylinders (gas flow velocity, 1 L/min). The C. putoensis seedlings in another climate chamber constituted the control (CK) group, which was quantitatively flushed with filtered air (without NO 2 ) at the same time. The chambers underwent a light/dark cycle with a photoperiod of 13 h and had a relative humidity of 60/50 ± 4% (day/night), with a temperature of 25/20 ± 3 • C (day/night). The control and NO 2 -treated seedlings (30 replicates in each treatment) were fumigated for 3 days (6 h per day), and then they recovered for 30 days [16].
The NO 2 concentration within the climate chamber containing leaves exposed to 1 L/min of air was measured with an NO 2 analyzer (model ML Series). After being treated with NO 2 , the seedlings were placed in an artificially controlled greenhouse under a natural simulation environment for 30 days of recovery. The environmental conditions of the greenhouse were as follows: room temperature, 25-28 • C; relative humidity, 60-70%; photoperiod, 14 h; and photosynthetically active radiation, 1000 µmol photons/(m 2 s).
For the following experiments, whole leaves were used unless otherwise specified.

Determination of Total Peroxidase (POD) Activity
POD is a class I oxidation-reduction enzyme that acts as a catalyst in a variety of biological processes; thus, it is an essential protective enzyme against reactive oxygen cell damage [18]. In response to adversity, POD is activated and provides resistance against adverse oxidation stress [19]. In this study, the POD level was measured with a guaiacol colorimeter [20]. The samples were pooled, and approximately 0.2 g of fresh leaves was placed in a pre-chilled mortar and then ground with 0.2 g of quartz sand. A total of 6 mL of 0.05 mol/L phosphate buffer (pH, 7.5) was added (in three applications, including one for mortar rinsing). The resulting homogenate was poured into a 10 mL centrifuge tube and stored at 4 • C. Centrifugation was performed at 5000× g for 20 min, and the obtained supernatant was a crude extract of POD. The reaction system for measuring the enzymatic activity contained 2.9 mL of phosphate buffer (0.05 mol/L), 1.0 mL of H 2 O 2 (2%), 1.0 mL of guaiacol (0.05 mol/L), and 0.1 mL of enzymatic solution. The enzymatic solution was boiled for 5 min and used as the control. After the enzymatic solution was applied, the system was immediately subjected to a 15-min incubation at 37 • C, which was followed by an ice bath. Trichloroacetic acid (20%, 2.0 mL) was added to terminate the reaction. Filtration (Steripak-GP, 10 L; Millipore, Germany) and appropriate dilution were then performed. The absorbance was measured at 470 nm [20]. Six replicates were designed for each group.

Transmission Electron Microscopy (TEM)
The plant material was cut into 1-mm 2 pieces and then fixed with 2.5% glutaraldehyde in a 0.1 M sodium cacodylate buffer (pH 7.4) for 4 h. After three washes with cacodylate buffer, the samples were fixed in 2% (w/v) osmium tetroxide in cacodylate buffer for 2 h. The samples were embedded in epoxy resin and dehydrated with an acetone series. Sections were cut using an LKB III ultramicrotome at 1 µm for light microscopy (LM) and approximately 50 nm for TEM. Ultrathin sections were stained with uranyl acetate and basic lead citrate and then analyzed by a Hitachi Hu 12a electron microscope [16].

RNA Isolation, cDNA Library Construction, and Illumina Sequencing
To understand the changes in gene levels after NO 2 fumigation, we selected the CK group and the 72-h NO 2 treatment group for transcriptome sequencing analysis. Two groups were prepared: a NO 2 treatment group and a CK group. After the leaves were removed from the tree, they were pooled and immediately frozen in liquid nitrogen and then stored at −80 • C in an ultra-low temperature freezer. The total RNA was extracted using the cetyltrimethy lammonium bromide (CTAB) method [21] and treated with RNasefree DNase I (TaKaRa, Dalian, China). The total RNA integrity was checked using gel electrophoresis, and the content was quantified using an ND-1000 spectrophotometer (Thermo, Waltham, MA, USA). Oligo (dT) 25 magnetic beads were used for isolating poly-(A) tail-containing mRNAs from the total RNA (20 µg), and mRNA was disrupted into short fragments with a fragmentation buffer at 70 • C for 5 min. These short fragments were used as templates to synthesize first-strand cDNA using random hexamer primers and reverse transcriptase. Second-strand cDNA fragments were obtained using a buffer containing DNA polymerase I, dNTPs, and RNase H. The final cDNA library was obtained by ligating the cDNA fragments to sequencing adaptors (Genomic DNA Sample Preparation Kit, Illumina, San Diego, CA, USA; two terminal sequencing: read length, 150 bp; paired end) and by conducting PCR amplification (Illumina Genomic Sample Preparation Kit, Illumina, San Diego, CA, USA). An Illumina HiSeq 2000 platform (Macrogen Bioinformatics Technology, Shenzhen, China) was used to sequence the mRNAs. Three replicates were designed for each group.

Data Analysis for RNA-seq Experiments
Adaptor sequences and low-quality reads were removed from the raw reads to obtain clean data [22,23]. The trinity method was adopted to assemble the clean data into transcripts [24]. National Center for Biotechnology Information, U.S. National Library of Medicine (NCBI) BLAST was used to compare the transcripts with NR, Swiss-Prot, Gene Ontology (GO), euKaryotic Orthologous Groups (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and several PFAM databases to obtain functional annotations [25]. The procedures for the RNAseq sequencing evaluation were as follows: Bowtie2 was used to compare the effective data from the samples to the spliced transcripts, and the mapping information was counted; Rseqc was used to analyze the redundant sequences and the distribution of inserted fragments; and BEDtools was used to check the homogeneity distribution and analyze the gene coverage [26]. A gene structure analysis was then performed. Specifically, BCFtools was used to seek possible SNP sites according to the mapping results; MISA was used for SSR analysis based on the sequence information of the spliced transcripts [27]. Salmon was used to calculate the gene expression. WGCNA was used for gene co-expression analysis. Based on the expression matrix of the samples, multi-directional statistical analyses and exploration, such as comparative analyses of the samples, were performed [28,29].

Identification, Annotation, and Enrichment Analysis of Differentially Expressed Genes
To identify differentially expressed genes (DEGs) related to the leaf metabolism of C. putoensis after NO 2 stress, we used RNA-seq by expectation maximization (RSEM) to map the clean reads of each sample to the transcriptome assemblies, and we used the DESeq with the following thresholds for DEG identification: false discovery rate (FDR), 0.01; fold change, 2 [30]. The identified DEGs were then used for GO and KOG classification and KEGG pathway enrichment analysis.

Validation by RT-qPCR
The results from the RNA-seq experiment were validated by analyzing eight plant genes that were most significantly differentially regulated under NO 2 stress (the smallest p-value was 1 × 10 −30 for chloroplasts) using RT-qPCR with cDNA as the template. RNA was obtained using the same method described in the Section 2.5. Oligo 7 software was used to design all the primers for RT-qPCR (Supplementary Table S1). A TB Green Premix Ex Taq kit (TaKaRa, Shiga, Japan) was used to perform RT-qPCR and an ABI StepOne plus thermal cycler (Applied Biosystems, Foster City, CA, USA) was used to run the RT-qPCR.

Morphology and Cell Structure of C. putoensis Leaves
The leaf morphology exhibited various changes when C. putoensis was exposed to NO 2 gas. According to Figure 1, the C. putoensis leaf damage appeared mostly as necrotic spots, from black spots to yellow spots, to an increasing extent. Some areas (such as the leaf tip) were severely damaged under NO 2 stress for 1-72 h. Figure  ences were observed between the plastids of the CK group and those of the NO2-treated plants that had recovered for 72 h (Figure 2 a,f).

Changes in POD Activity
Changes in the POD activity of C. putoensis at different NO 2 stress time points are shown in Figure 3. With increasing NO 2 fumigation time, the POD activity of C. putoensis increased, ranging from 385 U/(g min) fw to 596 U/(g min) fw. The 72-h treatment group had the highest POD level, with a significant difference compared to any of the remaining groups. Compared with the CK group, the 24 h treatment group showed a significant difference. The recovery group did not show a significant difference from the CK group.

Changes in POD Activity
Changes in the POD activity of C. putoensis at different NO2 stress time points are shown in Figure 3. With increasing NO2 fumigation time, the POD activity of C. putoensis increased, ranging from 385 U/(g min) fw to 596 U/(g min) fw. The 72-h treatment group had the highest POD level, with a significant difference compared to any of the remaining groups. Compared with the CK group, the 24 h treatment group showed a significant difference. The recovery group did not show a significant difference from the CK group.

RNA-seq Analysis of Clean Data from C. putoensis
C. putoensis is a non-model organism; therefore, de novo assembly is the only option for sequence assembly. In de novo assemblies, without the guidance of a reference sequence, the reads are assembled into contigs. To cover the C. putoensis transcripts completely, de novo assembly was used to generate the consensus transcriptome using Illumina sequencing data from samples under two different conditions together with raw reads from NO2-treated leaves and CK leaves. Due to trimming (extra bases whose lengths were lower than 20) and duplicate removal, we analyzed 529,540 transcripts with an average length of 425.97 bp for the de novo assembly of 250,200 unigenes with an average length of 376.73 bp (Table 1). In total, the highest annotation ratio was achieved for the GO database (110,530, 44.18%) ( Table 2), which represents successful annotation with known proteins. Only 1.84% of the genes were successfully annotated in all the databases; thus, many genes were without annotation. In this study, we focused on the sequence with the highest annotation ratio compared to the GO library to obtain the phase of the gene sequence and functional information for C. putoensis and its related species, as long as the gene over 136 K had at least 1 annotation. According to the GO classification (Figure 4a  3.3. RNA-seq Analysis of Clean Data from C. putoensis C. putoensis is a non-model organism; therefore, de novo assembly is the only option for sequence assembly. In de novo assemblies, without the guidance of a reference sequence, the reads are assembled into contigs. To cover the C. putoensis transcripts completely, de novo assembly was used to generate the consensus transcriptome using Illumina sequencing data from samples under two different conditions together with raw reads from NO 2 -treated leaves and CK leaves. Due to trimming (extra bases whose lengths were lower than 20) and duplicate removal, we analyzed 529,540 transcripts with an average length of 425.97 bp for the de novo assembly of 250,200 unigenes with an average length of 376.73 bp (Table 1). In total, the highest annotation ratio was achieved for the GO database (110,530, 44.18%) ( Table 2), which represents successful annotation with known proteins. Only 1.84% of the genes were successfully annotated in all the databases; thus, many genes were without annotation. In this study, we focused on the sequence with the highest annotation ratio compared to the GO library to obtain the phase of the gene sequence and functional information for C. putoensis and its related species, as long as the gene over 136 K had at least 1 annotation. According to the GO classification (Figure 4a), biological processes (274,614 genes, 36.98%), cellular components (236,419 genes, 31.84%), and molecular functions (231,488 genes, 31.176%) were identified. The KOG classification included 25 functional categories, including posttranslational modification, protein turnover, chaperones (7858 genes, 12.23%), translation, ribosomal structure and biogenesis (6309 genes, 9.82%), and general function prediction only (7041 genes, 10.96%) (Figure 4b). Additionally, the annotated genes were enriched in 23 KEGG pathways (Figure 4c). The top six enriched pathways included translation, carbohydrate metabolism, signal transduction, folding sorting and degradation, overview, and amino acid metabolism. tabolism, signal transduction, folding sorting and degradation, overview, and amino acid metabolism.

Identification and Analysis of DEGs in C. putoensis Leaves under NO 2 Stress
As in the experimental chambers, all the physical parameters other than the NO 2 concentration were kept the same; therefore, we presume that the observed results are solely caused by elevated NO 2 . Through the analysis of the CK group and the NO 2 stress group, the regulatory mechanisms and key genes of C. putoensis NO 2 stress were further explored. To identify DEGs between the two different samples, we analyzed the genes expressed in the two groups; a Venn diagram showed the distribution of specific genes (79,437 and 70,248 expressed genes in the control group (A) and the stressed group (B), respectively) and shared genes (99,724 expressed genes) ( Figure 5). Afterwards, pairwise comparisons are performed with FC ≥ 2 and FDR < 0.01 as the standards. In total, the RNAseq data involved one pairwise comparison, and 2,359 DEGs were ultimately identified, including 1,119 upregulated genes and 1240 downregulated genes ( Table 3). The DEGs were annotated using the KOG (877 DEGs, 37.18%), GO (1686 DEGs, 71.47%), KEGG (277 DEGs, 11.74%), and NR (1830 DEGs, 76.6%) databases and the conserved domains database (CDD, 2359 DEGs, 100%) ( Table 3). A pairwise comparison of the volcano plots map clearly shows the distribution of upregulated and downregulated genes (Figure 6a). Transcription factors (TFs) are the key components of regulatory systems that control and modulate stress adaptive pathways [22]. In accordance with the highly significant roles of TFs under NO 2 stress, we analyzed all the genes to identify the top 30 TF families (Figure 6b), which predominantly included C2H2, Zn-clus, C3H, bZIP, AP2/ERF-ERF, GRAS, bHLH, MYB-related, WRKY, and NAC.

Identification and Analysis of DEGs in C. putoensis Leaves under NO2 Stress
As in the experimental chambers, all the physical parameters other than the NO2 concentration were kept the same; therefore, we presume that the observed results are solely caused by elevated NO2. Through the analysis of the CK group and the NO2 stress group, the regulatory mechanisms and key genes of C. putoensis NO2 stress were further explored. To identify DEGs between the two different samples, we analyzed the genes expressed in the two groups; a Venn diagram showed the distribution of specific genes (79,437 and 70,248 expressed genes in the control group (A) and the stressed group (B), respectively) and shared genes (99,724 expressed genes) ( Figure 5). Afterwards, pairwise comparisons are performed with FC ≥ 2 and FDR < 0.01 as the standards. In total, the RNA-seq data involved one pairwise comparison, and 2,359 DEGs were ultimately identified, including 1,119 upregulated genes and 1240 downregulated genes ( Table 3). The DEGs were annotated using the KOG (877 DEGs, 37.18%), GO (1686 DEGs, 71.47%), KEGG (277 DEGs, 11.74%), and NR (1830 DEGs, 76.6%) databases and the conserved domains database (CDD, 2359 DEGs, 100%) ( Table 3). A pairwise comparison of the volcano plots map clearly shows the distribution of upregulated and downregulated genes (Figure 6a). Transcription factors (TFs) are the key components of regulatory systems that control and modulate stress adaptive pathways [22]. In accordance with the highly significant roles of TFs under NO2 stress, we analyzed all the genes to identify the top 30 TF families (Figure 6b), which predominantly included C2H2, Zn-clus, C3H, bZIP, AP2/ERF-ERF, GRAS, bHLH, MYB-related, WRKY, and NAC.     CDD  KOG  GO  KEGG  NR  NT  Upregulated genes  1119  1119  330  690  91  740  597  Downregulated genes  1240  1240  547  996  186  1090  760  Total  2359  2359  877  1686  277 1830 1357  The most common enriched pathways were found under GO classification, KEGG pathways, and KOG enrichment. In this study, we analyzed the GO classification of upregulated and downregulated annotated DEGs and selected the 30 with the smallest Q value for a scatter plot of pathway enrichment (Figure 7). The upregulated genes were assigned to 30 biological pathways functionally. The top three upregulated genes were involved in multicellular organism development (GO: 0007275), plastids (GO: 0009536), and chloroplasts (GO: 0009507), and the downregulated genes predominantly reflected response to stimulus (GO: 0050896), response to stress (GO: 0006950), and oxidoreductase activity (GO: 0016491). We also analyzed 91 upregulated and 187 downregulated KEGG pathways annotated with DEGs and chose the 30 with the smallest Q values for scatter plots of the pathway enrichment (Figures 8 and 9). The upregulated genes were functionally assigned to 76 biological pathways; the top upregulated genes were involved in photosynthesis (ko00195) (Figure 8), and the downregulated genes predominantly represented the biosynthesis of amino acids (ko01230) and carbon metabolism (ko01200) (Figure 9). The KEGG pathways showed that the DEGs of the NO 2 -treated group were significantly related to photosynthesis (Figure 10), i.e., four differentially expressed genes were involved in photosynthesis in C. putoensis under NO 2 stress. Combined with the genes classified by GO, which involved plastids and chloroplasts, this finding is consistent with the observed leaf changes in C. putoensis under NO 2 stress, i.e., the color change from green to yellow (shown in Figure 1). This result is also consistent with the change in cell ultrastructure as the chloroplast gradually deforms and more plastid granules appear with increasing NO 2 stress treatment time, which is a type of abiotic stress (Figure 2).

DEGs DEG Number
The most common enriched pathways were found under GO classification, KEGG pathways, and KOG enrichment. In this study, we analyzed the GO classification of upregulated and downregulated annotated DEGs and selected the 30 with the smallest Q value for a scatter plot of pathway enrichment (Figure 7). The upregulated genes were assigned to 30 biological pathways functionally. The top three upregulated genes were involved in multicellular organism development (GO: 0007275), plastids (GO: 0009536), and chloroplasts (GO: 0009507), and the downregulated genes predominantly reflected response to stimulus (GO: 0050896), response to stress (GO: 0006950), and oxidoreductase activity (GO: 0016491). We also analyzed 91 upregulated and 187 downregulated KEGG pathways annotated with DEGs and chose the 30 with the smallest Q values for scatter plots of the pathway enrichment (Figures 8 and 9). The upregulated genes were functionally assigned to 76 biological pathways; the top upregulated genes were involved in photosynthesis (ko00195) (Figure 8), and the downregulated genes predominantly represented the biosynthesis of amino acids (ko01230) and carbon metabolism (ko01200) (Figure 9). The KEGG pathways showed that the DEGs of the NO2-treated group were significantly related to photosynthesis (Figure 10), i.e., four differentially expressed genes were involved in photosynthesis in C. putoensis under NO2 stress. Combined with the genes classified by GO, which involved plastids and chloroplasts, this finding is consistent with the observed leaf changes in C. putoensis under NO2 stress, i.e., the color change from green to yellow (shown in Figure 1). This result is also consistent with the change in cell ultrastructure as the chloroplast gradually deforms and more plastid granules appear with increasing NO2 stress treatment time, which is a type of abiotic stress (Figure 2).

RT-qPCR Analysis of NO2 Stress-related Genes
To calculate the accuracy of the RNA-seq, we selected the DEGs with the most significant differences related to NO2 stress. We used a functional prediction of annotated genes from the RNA-seq data to identify eight DEGs, namely TRINITY_DN86073_c6_g3 (peroxidase 12-like, POD1), TRINITY_DN80077_c8_g2 (allene oxide synthase, HPL1), TRIN-ITY_DN80077_c8_g3 (allene oxide synthase, HPL2), TRINITY_DN86773_c3_g1 (allene oxide synthase, HPL3), TRINITY_DN81001_c0_g2 (hypothetical protein CICLE, APX5), TRINI-TY_DN86877_c1_g5 (geranylgeranyl diphosphate reductase, chloroplastic, CHL2), TRINI-TY_DN84191_c2_g1 (chloroplast chlorophyll a/b binding protein, CHL3), and TRINI-TY_DN86070_c0_g3 (hypothetical protein, CHLA) ( Figure 11). RT-qPCR analysis was performed on 8 candidate genes to verify the expression pattern of RNA-seq data ( Figure  12). The differential expression profiles of DEGs were consistent between the RNA-seq and RT-qPCR data, except for those of CHLA. Although there is a significant difference in the expression profile of one gene, when the RT-qPCR data are compared with the RNA-seq data, there are seven genes that show similar expression profiles. Our study found that CHL2, CHL3, and CHLA genes showed lower expression levels in the C. putoensis leaves upon NO2 stress. Strikingly, the selected oxidation family genes POD1, HPL1, and APX5 exhibited higher expression in C. putoensis upon NO2 stress. These findings seem to suggest that these genes participate in regulating the physiological response of C. putoensis.

RT-qPCR Analysis of NO 2 Stress-related Genes
To calculate the accuracy of the RNA-seq, we selected the DEGs with the most significant differences related to NO 2 stress. We used a functional prediction of annotated genes from the RNA-seq data to identify eight DEGs, namely TRINITY_DN86073_c6_g3 (peroxidase 12-like, POD1), TRINITY_DN80077_c8_g2 (allene oxide synthase, HPL1), TRIN-ITY_DN80077_c8_g3 (allene oxide synthase, HPL2), TRINITY_DN86773_c3_g1 (allene oxide synthase, HPL3), TRINITY_DN81001_c0_g2 (hypothetical protein CICLE, APX5), TRINITY_ DN86877_c1_g5 (geranylgeranyl diphosphate reductase, chloroplastic, CHL2), TRINITY_ DN84191_c2_g1 (chloroplast chlorophyll a/b binding protein, CHL3), and TRINITY_DN86070_c0 _g3 (hypothetical protein, CHLA) ( Figure 11). RT-qPCR analysis was performed on 8 candidate genes to verify the expression pattern of RNA-seq data ( Figure 12). The differential expression profiles of DEGs were consistent between the RNA-seq and RT-qPCR data, except for those of CHLA. Although there is a significant difference in the expression profile of one gene, when the RT-qPCR data are compared with the RNA-seq data, there are seven genes that show similar expression profiles. Our study found that CHL2, CHL3, and CHLA genes showed lower expression levels in the C. putoensis leaves upon NO 2 stress. Strikingly, the selected oxidation family genes POD1, HPL1, and APX5 exhibited higher expression in C. putoensis upon NO 2 stress. These findings seem to suggest that these genes participate in regulating the physiological response of C. putoensis. Figure 11. RT-qPCR validations of 8 candidate genes involved in NO2 stress in C. putoensis based on RNA-seq data. Hypothetical protein, chloroplastic, peroxidase, and allene oxide synthase represent different gene types.

Discussion
The results show that gaseous NO2 has a significant impact on the ultrastructure of mesophyll cells, i.e., increased translucence in the plastoglobuli, decreased chloroplasts and an increased number of plastoglobuli. Compared with that of the control group, the results are consistent with the gaseous SO2 and NO2 that cause swelling in the thylakoids and a decrease in the number of grana stacks [31]. The observed changes in the leaf cell structure are similar to those described in Ca-induced plants in the stressed group [30], namely, irregular plastid shape. Part of the reason for these changes may be that NO2 changes the semi-permeability of the plastid envelope. NO2 can interact directly with lipids, which is probably related to membrane effects [11]. The effects of chemical substances, such as H2O2 [32], ascorbic acid [33], and Na2S [34], have been studied before. However, the effects of natural restoration on plant responses to atmospheric pollution, especially NO2, has not been reported before. Our results indicate that natural recovery could be helpful for cell structure recovery and chloroplast morphology. No significant differences were observed between the CK group and the recovered plants, which is consistent with the findings of Souza et al. [35], who found that natural recovery from

Discussion
The results show that gaseous NO2 has a significant impact on the ultrastructure of mesophyll cells, i.e., increased translucence in the plastoglobuli, decreased chloroplasts and an increased number of plastoglobuli. Compared with that of the control group, the results are consistent with the gaseous SO2 and NO2 that cause swelling in the thylakoids and a decrease in the number of grana stacks [31]. The observed changes in the leaf cell structure are similar to those described in Ca-induced plants in the stressed group [30], namely, irregular plastid shape. Part of the reason for these changes may be that NO2 changes the semi-permeability of the plastid envelope. NO2 can interact directly with lipids, which is probably related to membrane effects [11]. The effects of chemical substances, such as H2O2 [32], ascorbic acid [33], and Na2S [34], have been studied before. However, the effects of natural restoration on plant responses to atmospheric pollution, especially NO2, has not been reported before. Our results indicate that natural recovery could be helpful for cell structure recovery and chloroplast morphology. No significant differences were observed between the CK group and the recovered plants, which is consistent with the findings of Souza et al. [35], who found that natural recovery from

Discussion
The results show that gaseous NO 2 has a significant impact on the ultrastructure of mesophyll cells, i.e., increased translucence in the plastoglobuli, decreased chloroplasts and an increased number of plastoglobuli. Compared with that of the control group, the results are consistent with the gaseous SO 2 and NO 2 that cause swelling in the thylakoids and a decrease in the number of grana stacks [31]. The observed changes in the leaf cell structure are similar to those described in Ca-induced plants in the stressed group [30], namely, irregular plastid shape. Part of the reason for these changes may be that NO 2 changes the semi-permeability of the plastid envelope. NO 2 can interact directly with lipids, which is probably related to membrane effects [11]. The effects of chemical substances, such as H 2 O 2 [32], ascorbic acid [33], and Na 2 S [34], have been studied before. However, the effects of natural restoration on plant responses to atmospheric pollution, especially NO 2 , has not been reported before. Our results indicate that natural recovery could be helpful for cell structure recovery and chloroplast morphology. No significant differences were observed between the CK group and the recovered plants, which is consistent with the findings of Souza et al. [35], who found that natural recovery from water stress could lead to the complete recovery of all gas exchange three days after rewatering.
As an important antioxidant enzyme, POD scavenges reactive oxygen species (ROS) [36]. In our experiment, the POD activity increased under NO 2 stress, indicating that C. putoensis plants exhibit substantial ROS-scavenging ability under NO 2 stress. In tolerant plant species, POD activity is higher, which enables the plants to protect themselves against oxidative stress [37,38]. In C. putoensis, it is not known how these changes at the cellular level are regulated at the genetic level. Therefore, we selected the CK group and 72-h NO 2 treatment group for transcriptome sequencing analysis.
According to the highly significant role of TFs under NO 2 stress, we analyzed all the genes to identify the top 30 TF families (Figure 6b), which predominantly included C2H2, Zn-clus, C3H, bZIP, AP2/ERF-ERF, GRAS, bHLH, MYB-related, WRKY, and NAC. These TF families are widely present in a variety of plant species, and they participate in the control of plant development and responses to biotics and biotic stress [39]. Previous research has revealed only the complete chloroplast genome of C. putoensis [40]. Our study is the first exploration of these TF families in C. putoensis based on transcriptome analysis.
In our experiment, many types of TFs, such as bZIP, NAC, AP2/ERF, and MYB, are involved in drought stress responses, and AP2/ERF-ERF is a large family of TFs in plants. AP2/ERF-ERF TFs are identified by the presence of an AP2 DNA-binding domain composed of 60-70 highly conserved amino acids. AP2/ERF-ERF TFs have significant functions in biological processes, including development, reproduction, primary and secondary metabolite biosynthesis, and adaptation to biotic and abiotic stresses [41]. They are primarily activated in response to drought stress [42], heat [43], waterlogging [44], high salinity [45], and osmotic stress [46]; however, this study is the first example of their activation in response to NO 2 stress. According to the literature, MYB TFs play a role in metabolism, cell fate and identity, development, and responses to biotic and abiotic stresses during the plant life cycle [47]. The roles of WRKY TFs in plant development, hormone signaling, biotic stress, and abiotic stress have been demonstrated [48]. A transcriptome analysis of Arabidopsis roots also indicated the upregulation and downregulation of WRKY TFs in response to NO 2 stress [49]. Plant-specific NAC transcription factors have multiple functions, including plant development, defense, and abiotic stress [50]; different plants have different abiotic stress responses to NAC-TFs [51]. However, all of the above TFs were determined to have roles in a NO 2 stress response, and NO 2 exposure is a type of abiotic stress.
In our study, several genes that were induced coded for photosynthesis-antenna proteins, and this expression was altered, as shown in Figure 11. The reduction in photosynthesis may also be attributed to degradation and damage in the thylakoid membrane protein-pigment complexes, and possibly also effects on lipids, thereby inducing oxidative stress in stressed plants [52]. As part of a defense mechanism to reduce the oxidative stressed damage, scavenging enzymes such as POD may be activated [53]. In our research, the expression levels of several genes for these enzymes and proteins were modulated. The role of these antioxidants includes altering gene expression to provide a redox buffer and act as a metabolic interface to regulate the optimum induction of adaptive responses [54]. NO 2 stress has adverse effects on plant growth and productivity; in higher plants, the photosynthesis apparatus is reorganized for acclimation to environmental and metabolic conditions [55]. However, reduced growth under stress is associated with an increase in photosynthesis-related genes, indicating sustained photosynthetic activity under NO 2 stress [56]. NO 2 stress leads to enhanced ROS production. In earlier reports, NO 2 treatment also significantly improved the antioxidant and isozyme activities, including those of superoxide dismutase and POD [57]. These enzymes catalyze the biosynthetic steps of various plant metabolites, and several researchers have demonstrated their role in stress tolerance [58]. Of the 87 POD genes, the majority were significantly upregulated after NO 2 stress, which is consistent with the increase in total POD activity [59]. This is a common response to various oxidative stress factors. Our study identified six differentially expressed transcripts encoding PODs, which are likely involved in the detoxification of ROS in C. putoensis under NO 2 stress and may be potential candidate genes for increasing NO 2 tolerance.
The results of our study indicate that the effects of NO 2 exposure on the ultrastructure of the cell structure, POD activity, and morphological changes were directly related to the NO 2 treatment time; therefore, we speculate that the effects of NO 2 on plants is partly attributable to the generation and accumulation of N 2 -derived NO 2 − in apoplastic and symplastic spaces. Moreover, with increasing NO 2 exposure, C. putoensis leaves were stressed by a large amount of NO 2 within a short time, followed by cell membrane damage and chloroplast destruction; this destruction then affected leaf photosynthesis, which altered the expression of genes related to abiotic stress.
Understanding plant stress responses to gaseous pollution is very important for urban greening applications. This study represents the first transcriptome analysis of NO 2 stress in C. putoensis, which is a relatively new field of research, particularly regarding photosynthesis and redox aspects. Therefore, RNA-seq analysis should urgently be conducted to provide further insights into these processes.

Conclusions
In this study, we recorded the changes in the morphology and anatomy of C. putoensis leaves under NO 2 stress. NO 2 stress adversely affected the morphology, leaf anatomy, and POD activity. These findings extend our understanding of plant stress responses; they also strongly indicate a need for further RNA-seq analysis. In this study, we used NCBI to compare the transcript with nine databases and obtained its functional annotations. We annotated 2,255 million clean Illumina paired-end RNA-seq reads (clean means to remove the bases with the mass value of reads less than 20 from raw data), and 250,200 unigene sequences were assembled based on the transcriptome data, with an average length of 376.73 bp and an N50 of 381 bp. A comprehensive functional annotation provided functional descriptions for more than 89% of the C. putoensis transcripts. Under NO 2 stress treatment, the plants had 2359 DEGs, among which 1119 exhibited upregulated expression and 1240 exhibited downregulated expression. GO enrichment analysis showed that the DEGs predominantly involved substance metabolism, protein binding, and catalytic activity. The DEGs were typically involved in metabolic pathways and photosynthesis metabolism, which were presented by KEGG analysis. According to the KOG analysis, DEGs were predominantly involved in carbohydrate transport and metabolism, translation, ribosome structure and biogenesis, biosynthesis, and the transport and catabolism of secondary metabolites. According to the KEGG pathway analysis, the expression of photosynthetic genes may be affected by NO 2 stress. Moreover, GO classification analysis indicated that the chloroplasts, plastids, and stimulus response may be related to NO 2 stress. Additionally, we found that the expression of POD families experienced dynamic changes during NO 2 stress treatment. According to the RT-qPCR validation, the HPL2, HPL1, and POD genes in C. putoensis also exhibited high expression under NO 2 stress. This study provides new insights into the C. putoensis processes that occur during NO 2 stress. Furthermore, the resulting transcriptome data represent an important candidate gene resource for future plant gene structure studies. These data will be very helpful during genome annotation in future genome projects.

Conflicts of Interest:
The authors declare no conflict of interest.