Abstract 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 NO[2] stress. The molecular mechanism of the C. putoensis NO[2] stress response was evaluated using sequencing data. NO[2] 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 NO[2] 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 NO[2] stress. Additionally, NO[2] 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 NO[2] stress supplies a reference for studying the molecular mechanism of C. putoensis resistance to NO[2] stress. The given transcriptome data represent a valuable resource for studies on plant genes, which will contribute towards genome annotations during future genome projects. Keywords: transcriptome, NO[2] stress, high-throughput sequencing, molecular mechanism, resistance, gene expression 1. 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 [[32]1]. At high temperatures, NO[2] is a maroon gas with a typically harsh odor, and it is a key contributor to air pollution [[33]2]. NO[2] is also an important component of acid rain [[34]3]. Its corrosivity and highly oxidative nature make it harmful to plant biochemical and physiological processes after entering plants through the stomata [[35]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 [[36]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) [[37]6], which can be incorporated into the α-amino group of soluble free amino acids [[38]7,[39]8], thereby not causing harm to the leaves [[40]9,[41]10]. Mayer et al. [[42]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 [[43]11]. The second theory is that the majority of plants have a low absorption capacity for NO[2]-N incorporation into the total plant N and can resist NO[2] [[44]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. 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 [[45]13]. The Zhejiang Forestry Science Research Institute has researched the cultivation and breeding of C. putoensis [[46]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 [[47]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 [[48]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. 2. Materials and Methods 2.1. 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. [[49]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 [[50]16]. 2.2. NO[2] Fumigation Fumigation was performed according to the method described in the literature [[51]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 [[52]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. 2.3. 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 [[53]18]. In response to adversity, POD is activated and provides resistance against adverse oxidation stress [[54]19]. In this study, the POD level was measured with a guaiacol colorimeter [[55]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 [[56]20]. Six replicates were designed for each group. 2.4. 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 [[57]16]. 2.5. 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 [[58]21] and treated with RNase-free 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. 2.6. Data Analysis for RNA-seq Experiments Adaptor sequences and low-quality reads were removed from the raw reads to obtain clean data [[59]22,[60]23]. The trinity method was adopted to assemble the clean data into transcripts [[61]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 [[62]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 [[63]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 [[64]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 [[65]28,[66]29]. 2.7. 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 [[67]30]. The identified DEGs were then used for GO and KOG classification and KEGG pathway enrichment analysis. 2.8. 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 [68]Section 2.5. Oligo 7 software was used to design all the primers for RT-qPCR ([69]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. 3. Results 3.1. 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 [70]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. [71]Figure 2 shows the micrographs of leaves from the control group and the treated plants after 1 h, 6 h, 24 h, and 72 h of NO[2] treatment under a TEM. The CK group exhibited oval cells with numerous well-compartmentalized grana stacks within [72]Figure 2a. Minimal starch was present. In NO[2]-treated plants ([73]Figure 2b,c), most cell structures were slightly damaged. Within some cells, cell dehydration and shrinkage ([74]Figure 2d) and chloroplast deformation ([75]Figure 2e) were found. Plastoglobuli were more numerous in the disrupted cells. The chloroplasts of the recovery plant leaves are discoidal and rich in starch grains ([76]Figure 2f). Following recovery from the NO[2] treatment, almost all the cell structures and chloroplasts in the plants seemed to recover their normal morphology. Starch was present in all cases. In some chloroplasts, the thylakoid system was almost unaffected. Few differences were observed between the plastids of the CK group and those of the NO[2]-treated plants that had recovered for 72 h ([77]Figure 2a,f). Figure 1. [78]Figure 1 [79]Open in a new tab Morphological changes in C. putoensis leaves receiving NO[2] treatment. CK: control group. NO[2] treated: 1 h, 6 h, 24 h, 72 h, and recovery for 30 days. Figure 2. [80]Figure 2 [81]Open in a new tab Images of cell structures from the primary leaf under a transmission electron microscope. (a), control; (b–e), NO[2]-treated plants; and (f), recovery plants. From 1 h to 72 h, the plastoglobuli in the cells gradually increase, and the chloroplasts gradually shrink, and they become slender and sticky. Slowly, the cell wall is separated. The red arrows in (e) indicate plasmolysis. V, vacuole; P, plastoglobuli; PM, plasma membrane; S, starch grain; CW, cell wall; and Chl, chloroplasts. 3.2. Changes in POD Activity Changes in the POD activity of C. putoensis at different NO[2] stress time points are shown in [82]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. Figure 3. [83]Figure 3 [84]Open in a new tab POD activity after different NO[2] stress times or 30 days of recovery. Six replicates for each group. 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 ([85]Table 1). Table 1. Length and internal length distribution of transcripts and unigenes. Type No. ≥500 bp ≥1000 bp N50 N90 Maximum Length Minimum Length Total Length Average Length CG Content Transcript 529,540 124,713 29,088 470 231 7490 201 225,567,341 425.97 40–50% Unigene 250,200 41,790 9609 381 221 7490 201 94,258,132 376.73 60–70% [86]Open in a new tab In total, the highest annotation ratio was achieved for the GO database (110,530, 44.18%) ([87]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 ([88]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%) ([89]Figure 4b). Additionally, the annotated genes were enriched in 23 KEGG pathways ([90]Figure 4c). The top six enriched pathways included translation, carbohydrate metabolism, signal transduction, folding sorting and degradation, overview, and amino acid metabolism. Table 2. Annotation of unigenes in each database. Database Number of Genes Percentage (%) CDD 79,760 31.88 KOG 64,226 25.67 NR 94,267 37.68 NT 77,874 31.12 PFAM 51,696 20.66 Swiss-prot 103,389 41.32 TrEMBL 93,882 37.52 GO 110,530 44.18 KEGG 9284 3.71 At least one database 136,276 54.47 All database 4595 1.84 [91]Open in a new tab Figure 4. [92]Figure 4 [93]Open in a new tab GO (a), euKaryotic Ortholog Group (KOG) (b) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (c) classification of all the identified genes. 3.4. 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) ([94]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 ([95]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%) ([96]Table 3). A pairwise comparison of the volcano plots map clearly shows the distribution of upregulated and downregulated genes ([97]Figure 6a). Transcription factors (TFs) are the key components of regulatory systems that control and modulate stress adaptive pathways [[98]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 ([99]Figure 6b), which predominantly included C2H2, Zn-clus, C3H, bZIP, AP2/ERF-ERF, GRAS, bHLH, MYB-related, WRKY, and NAC. Figure 5. [100]Figure 5 [101]Open in a new tab Venn diagram analysis of the expressed genes in two samples (A: Control, B: NO[2] stressed). Table 3. Annotation of A vs B DEGs in a pairwise comparison. DEGs DEG Number 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 [102]Open in a new tab Figure 6. [103]Figure 6 [104]Open in a new tab Volcano plots (a) of RNA-seq data for an A vs B pairwise comparison; (b) the top 30 differentially-regulated TF families were identified among all the genes. A: NO[2] control. B: NO[2] stressed. q-value, <0.05; fold change, >2. 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 ([105]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 ([106]Figure 8 and [107]Figure 9). The upregulated genes were functionally assigned to 76 biological pathways; the top upregulated genes were involved in photosynthesis (ko00195) ([108]Figure 8), and the downregulated genes predominantly represented the biosynthesis of amino acids (ko01230) and carbon metabolism (ko01200) ([109]Figure 9). The KEGG pathways showed that the DEGs of the NO[2]-treated group were significantly related to photosynthesis ([110]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 [111]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 ([112]Figure 2). Figure 7. [113]Figure 7 [114]Open in a new tab GO enrichment factor analysis of the DEGs. (a) Upregulated genes: the top three upregulated genes are involved in multicellular organism development, plastids, and chloroplasts; (b) downregulated genes: the downregulated genes predominantly reflected response to stimulus, response to stress, and oxidoreductase activity. Figure 8. [115]Figure 8 [116]Open in a new tab Upregulated genes of KEGG pathway categories (a) and enrichment factor analysis (b) of the DEGs. The upregulated genes are functionally assigned to 76 biological pathways, and the top upregulated genes are involved in photosynthesis. Figure 9. [117]Figure 9 [118]Open in a new tab Downregulated genes of KEGG pathway categories (a) and enrichment factor analysis (b) of the DEGs. The downregulated genes predominantly represent amino acid biosynthesis and carbon metabolism. Figure 10. [119]Figure 10 [120]Open in a new tab There are 4 genes (Psb, Psa, Pet, and F-type ATPase a) involved in photosynthesis in C. putoensis under NO[2] stress. Green represents the downregulated expression of the gene, red represents the upregulated expression of the gene, and yellow indicates no significant difference in gene expression. 3.5. 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), TRINITY_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) ([121]Figure 11). RT-qPCR analysis was performed on 8 candidate genes to verify the expression pattern of RNA-seq data ([122]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. [123]Figure 11 [124]Open in a new tab RT-qPCR validations of 8 candidate genes involved in NO[2] stress in C. putoensis based on RNA-seq data. Hypothetical protein, chloroplastic, peroxidase, and allene oxide synthase represent different gene types. Figure 12. [125]Figure 12 [126]Open in a new tab The expression profiles according to RT-qPCR (relative expression) and RNA-seq (FPKM values). 4. 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 [[127]31]. The observed changes in the leaf cell structure are similar to those described in Ca-induced plants in the stressed group [[128]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 [[129]11]. The effects of chemical substances, such as H[2]O[2] [[130]32], ascorbic acid [[131]33], and Na[2]S [[132]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. [[133]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) [[134]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 [[135]37,[136]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 ([137]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 [[138]39]. Previous research has revealed only the complete chloroplast genome of C. putoensis [[139]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 [[140]41]. They are primarily activated in response to drought stress [[141]42], heat [[142]43], waterlogging [[143]44], high salinity [[144]45], and osmotic stress [[145]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 [[146]47]. The roles of WRKY TFs in plant development, hormone signaling, biotic stress, and abiotic stress have been demonstrated [[147]48]. A transcriptome analysis of Arabidopsis roots also indicated the upregulation and downregulation of WRKY TFs in response to NO[2] stress [[148]49]. Plant-specific NAC transcription factors have multiple functions, including plant development, defense, and abiotic stress [[149]50]; different plants have different abiotic stress responses to NAC-TFs [[150]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 [151]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 [[152]52]. As part of a defense mechanism to reduce the oxidative stressed damage, scavenging enzymes such as POD may be activated [[153]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 [[154]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 [[155]55]. However, reduced growth under stress is associated with an increase in photosynthesis-related genes, indicating sustained photosynthetic activity under NO[2] stress [[156]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 [[157]57]. These enzymes catalyze the biosynthetic steps of various plant metabolites, and several researchers have demonstrated their role in stress tolerance [[158]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 [[159]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. 5. 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. Supplementary Materials The following are available online at [160]https://www.mdpi.com/article/10.3390/genes12050754/s1, Table S1: List of 8 DEG primers used for RT-qPCR. Raw sequence data are available at [161]http://www.ncbi.nlm.nih.gov/bioproject/717124 (accessed on 10 May 2021), with BioSample accession number of SAMN19071592, SAMN19071593 and SRA accession number of PRJNA717124. [162]Click here for additional data file.^ (41.9KB, zip) Author Contributions Conceptualization, Q.S. and Z.Z.; methodology, Z.Z.; software, C.L.; validation, C.L. and M.S.; formal analysis, M.S.; investigation, Q.S. and C.L.; resources, Z.Z.; data curation, J.X.; writing, original draft preparation, C.L.; writing, review, and editing, Z.Z.; visualization, J.X.; supervision, Z.Z.; project administration, Z.Z.; and funding acquisition, Z.Z. All authors have read and agreed to the published version of the manuscript. Funding This study was supported by the National Natural Science Foundation of China (31770752), the 333 Projects of Jiangsu Province (BRA2018065), the Science and Technology Support Plan Project of Jiangsu Province (BM2013478), the Key Project of 2019 Teaching Quality Improvement Project: Research on the flipped classroom teaching mode of garden plant application, Jiangsu Forestry University, and the High Education Research Topics 2019: Research on the classroom reform mode of cultivating innovative landscape talents under the idea of ecological civilization (Jiangsu Forestry University). Institutional Review Board Statement Not applicable. Informed Consent Statement Not applicable. Data Availability Statement The data used for the analysis in this study are available within the article and [163]supplementary table. Conflicts of Interest The authors declare no conflict of interest. Footnotes Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. References