Ozone and nitrogen dioxide regulate similar gene expression responses in Arabidopsis but natural variation in the extent of cell death is likely controlled by different genetic loci

High doses of ozone (O3) and nitrogen dioxide (NO2) cause damage and cell death in plants. These two gases are among the most harmful air pollutants for ecosystems and therefore it is important to understand how plant resistance or sensitivity to these gases work at the molecular level and its genetic control. We compared transcriptome data from O3 and NO2 fumigations to other cell death related treatments, as well as individual marker gene transcript level in different Arabidopsis thaliana accessions. Our analysis revealed that O3 and NO2 trigger very similar gene expression responses that include genes involved in pathogen resistance, cell death and ethylene signaling. However, we also identified exceptions, for example RBOHF encoding a reactive oxygen species producing RESPIRATORY BURST OXIDASE PROTEIN F. This gene had increased transcript levels by O3 but decreased transcript levels by NO2, showing that plants can identify each of the gases separately and activate distinct signaling pathways. To understand the genetics, we conducted a genome wide association study (GWAS) on O3 and NO2 tolerance of natural Arabidopsis accessions. Sensitivity to both gases seem to be controlled by several independent small effect loci and we did not find an overlap in the significantly associated regions. Further characterization of the GWAS candidate loci identified new regulators of O3 and NO2 induced cell death including ABH1, a protein that functions in abscisic acid signaling, mRNA splicing and miRNA processing. The GWAS results will facilitate further characterization of the control of programmed cell death and differences between oxidative and nitrosative stress in plants.


Introduction
Ozone (O 3 ) and nitrogen dioxide (NO 2 ) are common air pollutants that occur in the troposphere, the lowest part of the Earth's atmosphere. In high concentrations, they cause damage to plants and animals (Middleton, 1961;Brunekreef and Holgate, 2002;Ainsworth, 2017). O 3 and nitrogen oxides are amongst the three most harmful air pollutants in terms of damage to ecosystems (The European Environment Agency, 2020). Low to intermediate ppb (parts per billion) levels of O 3 are harmful to plants (McGrath et al., 2015). However, NO 2 is less toxic than O 3 , causing reduced plant growth and lesion formation at concentrations in the high ppb to low ppm (parts per million) range (Taylor and Eaton, 1966;Kasten et al., 2016).
Several studies indicate large-scale yield loss due to O 3 pollution in agriculturally important species, including maize, soybean, rice and wheat (McGrath et al., 2015;Feng et al., 2022b).
The mechanisms leading to damage and cell death after exposure to high levels of O 3 are better understood than those of NO 2 . Both gases enter the plants through the stomatal pores. In the apoplastic space, O 3 forms reactive oxygen species (ROS), including superoxide and hydrogen peroxide (Vainonen and Kangasjarvi, 2015). Whether plants tolerate the amount of ROS produced by O 3 or initiate cell death depends on antioxidant capacity and the balance of stress hormone signaling. Increased production of ethylene promotes O 3 induced cell death in several plant species (Tuomainen et al., 1997;Overmyer et al., 2000;Vahala et al., 2003), whereas jasmonic acid (JA) protects from O 3 damage (Rao et al., 2000;Xu et al., 2015a). Salicylic acid (SA) is involved in both promotion of O 3 induced cell death as well as activation of defence response to O 3 (Rao and Davis, 1999;Xu et al., 2015b). The role of hormones in O 3 responses has been studied with A. thaliana mutants such as the JA receptor mutant coi1, and the ethylene overproducing mutant eto1, that are both O 3 sensitive (Rao et al., 2002;Xu et al., 2015a). After entering the apoplast, NO 2 rapidly reacts with water to produce nitrate, nitrite, nitric oxide (NO), and protons. Especially nitrite and NO are important drivers of NO 2 induced cell death in A. thaliana (Kasten et al., 2016). Mutant analyses suggested that signaling by SA enhanced NO 2 tolerance whereas JA does not regulate NO 2 induced cell death (Kasten et al., 2016). Small molecule antioxidants such as ascorbic acid (vitamin C) are important scavengers of ROS and RNS in plant tissues. Accordingly, the ascorbic acid deficient mutant vtc1 is sensitive to both O 3 as well as NO 2 (Conklin et al., 1996;Kasten et al., 2016).
Plants produce ROS and reactive nitrogen species (RNS) as second messengers to regulate growth and development, stress responses and long-distance signaling (Waszczak et al., 2018;Hancock and Neill, 2019;Castro et al., 2021). In plant-pathogen signaling, both ROS and NO acts as signal molecules, and are proposed to participate in an amplification loop to enhance signaling (Wang et al., 2018). Further, the balance between ROS and NO signaling regulates programmed cell death (PCD) (Delledonne et al., 2001;Wendehenne et al., 2014) Therefore, at least some effects of O 3 and NO 2 (or their reactive derivatives) relates to activation of cellular signaling processes. This assumption is supported by increased transcript levels of pathogen responsive and PCD related genes by both O 3 and NO 2 (Xu et al., 2015a;Mayer et al., 2018). Consequently, exposure to these gasses leads to the establishment of basal pathogen resistance in A. thaliana against the bacterial pathogen Pseudomonas syringae (Sharma et al., 1996;Mayer et al., 2018). Infection with avirulent pathogens elicit the hypersensitive defence response (HR) culminating in localized PCD, thought to restrict pathogen growth (Delledonne et al., 2001;Gaupels et al., 2011;Wang et al., 2013). O 3 -and NO 2 -triggered cell death is remarkably similar to HR-PCD since it is also dependent on simultaneous signaling by NO and ROS and is accompanied by the accumulation of fluorescent compounds in dying leaf tissues (Overmyer et al., 2005;Ahlfors et al., 2009;Kasten et al., 2016). Hence, O 3 and NO 2 cause leaf damage at least partially by activation of PCD (Gandin et al., 2021;Hussain et al., 2022).
Although the damaging effects of O 3 and NO 2 have been assessed in several studies, much less is known about the genetic factors controlling tolerance of plants to these pollutants. Regulation of O 3 and NO 2 responses in A. thaliana has mostly been studied in the genetic background Col-0 (Overmyer et al., 2008;Frank et al., 2019). Natural accessions of A. thaliana show large variation in their O 3 sensitivity . For instance, Cvi-0 is very sensitive but Col-0 rather tolerant to 300-350 ppb O 3 for 6 h . By studying naturally occurring variation for stress tolerance in multiple accessions, instead of a mutant approach in a single genetic background, it is possible to gain broader understanding of the underlying genetics. Natural variation can be explored with a genome wide association study (GWAS). This approach takes advantage of naturally occurring genetic recombination events to associate phenotypes of accessions with single nucleotide polymorphisms (SNPs) in close genomic vicinity of causative genes. The 1001 Genomes project has provided genome sequences of more than a thousand natural A. thaliana accessions and this data is available for use as genetic markers for GWAS (Alonso-Blanco et al., 2016).
To further understand how O 3 and NO 2 regulates defence signaling and cell death, we used two complementary approaches; transcriptome analysis and GWAS. The transcriptional responses to O 3 and NO 2 treatments were previously analysed individually (Brosche et al., 2014;Xu et al., 2015a;Mayer et al., 2018), and here we systematically identify the regulatory context of O 3 and NO 2 . In a comparison of different transcriptome datasets, we found very high overlaps of differentially expressed genes regulated by O 3 and NO 2 . With real time quantitative PCR (qPCR) we confirmed that both gasses activated similar signaling in different A. thaliana accessions. However, we also identified a marker gene with differential O 3 versus NO 2 response, demonstrating that A. thaliana can activate precise signal activation to each gas. We continued with GWAS for O 3 and NO 2 leaf damage with up to 372 A. thaliana natural accessions. We identified 12 genomic loci associated with O 3 and NO 2 induced leaf damage. Experiments with T-DNA knock-out mutants suggest functions of several GWAS-derived candidate genes in O 3 and NO 2 sensitivity.

Growth conditions
Experiments were performed in Helsinki and Munich. Plants in Helsinki were grown in 1:1 peat-vermiculite under 280 µmol m 2 s -1 white light irradiance, 12 h : 12 h light-darkcycle, 23°C/19°C (day/night) temperature and 70%/90% relative humidity. In Munich plants were grown in 5:1 Floradur propagation substrate -quartz sand, in walk-in size chambers, under 250 mmol m 2 s -1 ; PAR, 12 h : 12 h light-dark -cycle, 23°C/ 18°C (day/night) temperature and 70%/90% relative humidity. Plants grew faster in Munich than in Helsinki, possibly due to small differences in light and soil quality. However, plants were exposed to O 3 in Munich and Helsinki at a similar developmental stage (estimated by counting leaf numbers of the control plants Col-0 and Cvi-0). Therefore, plants were treated with O 3 approximately 4 days younger in Munich than in Helsinki.

O 3 fumigations in Helsinki
372 accessions at the age of 23 days were exposed to 400 ppb of O 3 for 6 has described previously . All treatments were started in the morning after approximately 2 h of light exposure. There were 8 plants of each genotype in one replicate. 23 accessions with controls could be treated simultaneously. Each genotype was present in two to three replicates. The O 3 phenotype was scored as number of injured leaves from all leaves relative to damage in Col-0.

O 3 fumigations in Munich
All 127 accessions were treated simultaneously in a single fumigation chamber. Five 19 days-old plants of every accession were fumigated for 6 h with 350 ppb of O 3 . The experiment was replicated twice.

NO 2 fumigations in Munich
Three 24-28 days-old plants of 216 accessions were exposed to NO 2 for 1 h in an air-tight fumigation chamber as reported earlier (Frank et al., 2019). Plants were short-term fumigated for 1 h with 10, 20 and 30 ppm of NO 2 . This range of concentrations was chosen as sensitive plants were damaged already with 10 ppm whereas tolerant plants displayed only weak symptoms even after 30 ppm of NO 2 . The NO 2 phenotype was scored as percent leaf area damaged (0% = score 1, <25% = score 2, 25-50% = 3, 50-75% = 4, >75% = 5). A cumulative score (scale from 3 to 15) from all three fumigations was used in further analyses. The experiment was repeated three times.

Analysis of public gene expression data
Hierarchical clustering was done with publicly available data (Supplementary Table S2). The raw data were processed with robust multiarray average normalization using Bioconductor limma and affy packages in R. Gene expression was summarized by calculating log 2 ratio of the treatment/mutant and control/wild type expression. Bayesian hierarchical clustering method was used (as described in (Wrzaczek et al., 2010)) with 1000 bootstrap resampling.
To compare similarities in O 3 and NO 2 transcriptomes we used data from exposures of A. thaliana Col-0 to 350 ppb O 3 for 2 h (RNA-seq analysis) (Xu et al., 2015a;Xu et al., 2015b) and 10 ppm NO 2 for 1 h (microarray analysis) (Mayer et al., 2018), or a NO experiment with the donor S-nitrosocysteine (1 mM, 6 h timepoint) (Hussain et al., 2016). Analysis of RNA-seq data is described in the Supplementary Methods. GO-term enrichment was performed in R (R Core Team 2018), version 3.5.0, using the package clusterProfiler (Yu et al., 2012). Venn diagrams were constructed using jvenn software (Bardou et al., 2014).

Real time reverse transcriptase quantitative PCR (qPCR)
Nine accessions were grown for qPCR analysis in Munich facilities in similar environmental and soil conditions as for GWAS. Plants were treated with 350 ppb O 3 or 10 ppm NO 2 at the age of three weeks. Three plants per genotype were collected and pooled in liquid nitrogen 2 h after the start of fumigations. RNA isolation, cDNA synthesis and qPCR were performed as described (Xu et al., 2015a). Normalization of the data was performed in qBase 2.0 [Biogazelle, (Hellemans et al., 2007)], with three reference genes SAND, TIP41 and YLS8 (Czechowski et al., 2005) The whole experiment was replicated three times. Primer sequences and amplification efficiencies can be found in Supplementary Table S3.

Genome-wide association studies
The GWA analyses were performed for the maximum numbers of phenotyped accessions for each trait, as well as separately for the set of accessions that was shared between all three experiments (119 accessions). The phenotype data from NO 2 treatment was nearly normally distributed but the phenotype data from both O 3 experiments were skewed as there were more tolerant accessions compared to sensitive. No transformations were used as these did not provide normality for the phenotypic data. Before the full genome data was available, the O 3 datasets were analysed with 250K SNP array data (Kim et al., 2007;Atwell et al., 2010) with EMMAX (Kang et al., 2010), from where some candidate genes were identified based on their biological function. The main GWA analyses of the damage screens were conducted in GWA-portal (http://gwas.gmi.oeaw. ac.at/, predecessor GWAPP by (Seren et al., 2012)), where imputed full genome data (Cao et al., 2011;Gan et al., 2011;Long et al., 2013), could be used for association analysis. The numbers of investigated SNPs were 4.1 million for 119 common accessions (O 3 and NO 2 ), 5.5 million for 372 accessions (O 3 Helsinki), 4.3 million for 127 accessions (O 3 Munich) and 4.9 million for 216 accessions (NO 2 ). The average SNP density in our GWA analyses was around one SNP per 25 base pairs, which should provide a very good coverage of the genome, even though linkage disequilibrium decays rapidly (on average within 10 kb (Kim et al., 2007)). We did not perform filtering based on minor allele frequency, but included all the SNPs in the analysis, as we wanted to include rare variants that play a role in O 3 sensitivity (Jakobson et al., 2016). Both non-parametric Kruskal-Wallis and accelerated mixed model (AMM) approaches were used. As variants truly associated with the traits of interest may occur in certain populations and therefore be correlated with population structure, we present data from Kruskal-Wallis tests, which does not correct for the population structure. AMM uses a linear mixed model approach that was developed by Kang et al. which corrects for population structure and genetic relatedness in association mapping (Kang et al., 2008;Kang et al., 2010).

Ion leakage of T-DNA mutant lines
To verify genomic regions identified by GWA, T-DNA lines (in Col-0 accession) of candidate genes were tested for O 3 and NO 2 sensitivity. As candidate genes we included all the genes that had SNPs in high linkage disequilibrium (r 2 > 0.8) with the Bonferroni corrected significant SNPs, as well as two genes that were selected based on their biological functions from the analysis with 250K SNP array data ( Table 1). The mutant lines are described in more detail in Supplementary Table S4. The T-DNA mutant lines were selected to have inserts in exons and were confirmed by PCR to be homozygous for the insert (PCR primers used for genotyping in Supplementary Table S4). Plants for ion leakage experiments were grown as described earlier and fumigated with O 3 (350 ppb, 6 h) in Helsinki and NO 2 (10 ppm for 6 h or 30 ppm for 1 h) in Munich. Ion leakage was performed as earlier (Brosche et al., 2014;Kasten et al., 2016). The experiments were repeated at least three times. Statistical analysis of the ion leakage measurements was done with linear mixed models in R 3.4.3 (R Core Team, 2018), with lme4 package. As several mutant lines were compared to the same control (Col-0), we used Dunnett's test, with multcomp package, to evaluate which comparisons were significant.

Results
3.1 O 3 and NO 2 trigger the expression of genes related to cell death and pathogen resistance ROS and RNS act as signals in plant stress responses, pathogen resistance, and PCD (Delledonne et al., 2001;Wang et al., 2013;Xu et al., 2015a;Mayer et al., 2018). To place O 3 -and NO 2 -induced transcriptome changes into the context of PCD, we performed Bayesian hierarchical clustering with all genes from the gene ontology (GO) category cell death ( Figure 1A). Transcriptome datasets included both air pollutants, cell death in lesion mimic mutants, pathogen infection and hormone treatments (for a full list of experiments see Supplementary  Table S2). The three resulting clusters ( Figure 1A) contained genes with strongly increased (cluster I), decreased (II), or weakly increased (III) transcript levels. In cluster I, there were large similarities between exposure to 350 ppb O 3 (2 h), 10 ppm NO 2 (1 h), pathogen infections (Botrytis cinerea, Pseudomonas syringae pv. maculicola) and lesion mimic mutants that undergo spontaneous cell death (acd11, mkk1 mkk2, siz1, ssi2). This highlights that the O 3 and NO 2 treatments trigger similar signaling pathways as those activated during pathogen infection and cell death.
To explore the overlap in transcriptional regulation between O 3 and NO 2 , we used transcriptome data from exposures of A. thaliana Col-0 to 350 ppb O 3 for 2 h (Xu et al., 2015b) and 10 ppm NO 2 for 1 h (Mayer et al., 2018) (Supplementary Table S5). In this comparison 1534 genes were >2-fold up-and 1016 genes down-regulated by both gases ( Figure 1B), i.e >55% of the genes regulated by NO 2 were also responsive to O 3 and vice versa. GO term enrichment analysis revealed that O 3 and NO 2 activated similar biological processes ( Figure 1C) including "regulation of immune response", "regulation of plant-type hypersensitive response", "respiratory burst", and "response to ethylene". In sum, O 3 and NO 2 transcriptionally regulate largely overlapping sets of genes involved in pathogen resistance, cell death, and ethylene signaling.
As NO is a main component in RNS signaling (Delledonne et al., 2001;Kasten et al., 2016), we also made transcriptome comparisons between O 3 and treatment with the NO-donor Snitrosocysteine (Hussain et al., 2016). Similar to the O 3 -NO 2 comparison, O 3 and S-nitrosocysteine regulated genes showed a large overlap with more than 6000 genes regulated by both treatments (Supplementary Figure S1A). GO enrichment showed that many aspects of defense signaling (including "phosphorelay signal transduction system", "response to ROS" and "regulation of response to biotic stimulus") were found among the genes with increased expression by both treatments (Supplementary Table S6).

O 3 -and NO 2 -induced defence transcriptional responses is similar between A. thaliana natural accessions, with notable exceptions
One drawback in the re-analysis of transcriptome data from different laboratories is that the experiments are heterogeneous in terms of plant growth conditions and treatments. Moreover, the transcriptome data shown in Figure 1 mainly originate from experiments with A. thaliana Col-0, and thus, might not reflect the natural variation within A. thaliana accessions. To address these issues, we grew and treated a set of accessions under controlled conditions followed by transcript level analysis using real time quantitative PCR (qPCR). The selected accessions covered different O 3 sensitivities (e.g. Ts-1 is O 3 tolerant and Cvi-0 O 3 sensitive), natural habitats, and genetic distances (Cvi-0 is more distantly related to other accessions [Alonso-Blanco et al., 2016)]. The plants were treated with 350 ppb O 3 and 10 ppm NO 2 as these treatments induced largely overlapping sets of genes in Col-0 ( Figure 1).
For qPCR analysis (Figure 2), we selected marker genes from the GO category cell death that also has been used as marker genes for the stress hormones SA, JA and ethylene. The ethylene/ JA marker genes CEJ1 (COOPERATIVELY REGULATED BY  Transcript levels after NO 2 and O 3 treatments in different accessions. Transcript levels of five marker genes was measured with qPCR in natural accessions of A. thaliana 2 h after treatments with NO 2 (left column) and O 3 (right column). Log 2 -transformed fold change of marker genes in natural accessions are shown from three biological replicates. Statistical analysis was performed with a linear mixed model (P < 0.05). Letters from a to e represent comparison of NO 2 treatment effect, and letters from A to D O 3 treatment effect, on different genotypes. Asterisks (*) represent statistical differences between NO 2 and O 3 treatments on each genotype. Error bars display standard deviation. The accessions were selected to include genotypes with differential NO 2 and O 3 damage phenotypes (see also Figure 3A). ETHYLENE AND JASMONATE 1) and RAP2.6 (RELATED TO AP2.6 (Krishnaswamy et al., 2011)) had increased transcript levels and were regulated similarly by both gases in all accessions. The SA marker gene GRX480 (GLUTAREDOXIN 480; (Caarls et al., 2015) and systemic signaling FMO1 (FLAVIN-DEPENDENT MONOOXYGENASE 1 (Hartmann et al., 2018), had increased transcript levels by both treatments, which were significantly higher by O 3 compared to NO 2 . Transcript levels for RBOHF (RESPIRATORY BURST OXIDASE HOMOLOG F), encoding a ROS producing enzyme, increased after O 3 but decreased after NO 2 treatment ( Figure 2). This represented the most contrasting effect of the two gases.
Overall, four of five studied marker genes had increased transcript levels by both pollutants with O 3 having a slightly stronger impact than NO 2 . There were no major differences between tolerant or sensitive accessions. For instance, the tolerant accession Ost-0 showed very similar transcript profiles as the sensitive accession Cvi-0. Taken together, Figures 1, 2 support the conclusion that O 3 and NO 2 trigger the expression of similar sets of defence-and cell death-related genes also in different A. thaliana accessions. However, the very contrasting regulation of RBOHF (increased transcript levels by O 3 and decreased transcript levels by NO 2 ), show that in addition to common signaling pathways induced by both gases, there is a mechanism by which A. thaliana can perceive and initiate signaling that is unique for O 3 versus NO 2 .

O 3 -and NO 2 -induced leaf symptoms are similar between accessions
Natural accessions show a wide range of O 3 tolerance versus sensitivity . As the transcriptome data ( Figure 1B, Supplementary Tables S5 and S6) has high overlap in O 3 and NO 2 induced biological mechanisms, we evaluated whether different accessions developed similar damage to O 3 and NO 2 . Ts-1, Ost-0, and Col-0 were tolerant while Bsch-0, Li-3, and Cvi-0 were sensitive to both pollutants. NFA-8 and Shahdara had intermediate phenotypes, e.g. Shahdara was clearly sensitive to O 3 in Helsinki conditions, less so in Munich and neither sensitive nor tolerant for NO 2 . Je-0 had a differential phenotype since it was rather tolerant to NO 2 but clearly sensitive to O 3 ( Figure 3A).
We extended the investigation of O 3 -and NO 2 -related leaf phenotypes to large sets of accessions (Supplementary Table S1). Independent O 3 damage screens took place in Helsinki where 372 accessions were treated for 6 h with 400 ppb O 3 and in Munich where 127 accessions received 350 ppb O 3 for 6 h. In both experiments results of the damage score (% of leaf area damaged) did not show normal distribution because tolerant accessions were overrepresented. In Munich, 216 accessions were treated with 10, 20, or 30 ppm NO 2 for 1 h. This treatment scheme facilitated a fine-tuned and nearly normally distributed rating of NO 2 leaf damage.
A subset of 119 accessions was common for all three experiments ( Figure 3B), which allowed us to compare damage from O 3 and NO 2 , as well as the comparison of O 3 damage at two different facilities. Altogether, the tested accessions exhibited large phenotypic variation after fumigation with O 3 or NO 2 (Figure 3, Supplementary Table S7). Scatter plots revealed a significant correlation in the extent of lesion formation between both O 3 experiments (Spearman's correlation coefficient r s = 0.45, P = 2.4 x 10 -7 , Figure 3C). Likewise, the phenotypes caused by the O 3 and NO 2 fumigations in Munich showed a good correlation (r s =0.32, P = 3.7 x 10 -4 , Figure 3D) although several accessions had differential sensitivities to both pollutants including Je-0 ( Figure 3A, see Supplementary Table S7 for further examples).

Identification of candidate genes involved in phenotypic responses to O 3 and NO 2 using GWAS
We used the leaf damage scores for the 119 accessions common to all three fumigation screens to perform GWAS, and to identify and compare SNPs associated with O 3 -and NO 2induced phenotypes. From GWAS, the P-values are a measure of association strength between SNPs and leaf damage of the different accessions. The analysis included both the nonparametric Kruskal-Wallis (KW) test and the accelerated mixed model (AMM) analysis that corrects for population structure. As a final outcome, this led to the identification of genes containing SNPs -or being within high linkage disequilibrium (LD) of SNPs -that were significantly associated with the damage score ( Figure 4 and Table 1). Overall, Bonferroni-adjusted AMM statistics identified three significant SNPs (with -log 10 P-values higher than 8) from the three O 3 and NO 2 experiments with the common 119 accessions ( Figure 4A). In the Helsinki dataset, a single associated SNP was found in the 5' UTR of AT3G61410 coding for a U-box kinase family protein. Two of these SNPs were significantly associated with O 3 induced leaf damage in the Munich dataset ( Figure 4A), where the smallest P-value SNPs were in an intergenic region closest to a gene with unknown function AT2G43795 and adjacent to MPK6 (MAP KINASE 6).
Subsequently, GWAS was performed on each individual dataset, with full number of accessions phenotyped, which identified additional SNPs significantly associated with O 3induced leaf phenotypes (Table 1) Supplementary Table S7 for the scores for each accession). O 3 sensitivity was quantified as O 3 -induced visible leaf injury (number of leaves with damage/total leaves, displayed as percentages and normalized to Col-0). The NO 2 phenotypes are visible injury scores from three fumigations (10, 20 and 30 ppm) with 3 plants per accession. The cumulative score is on a scale from 3 to 15. Spearman's correlation coefficients (r s ) are presented on the bottom right corner (C, D). The accessions Col-0 and Cvi-0, which were included in all treatments and repeats, are highlighted in blue and red, respectively. and one significant SNP upstream of AT1G44890. Localized between the latter two genes, MCM2 (MINICHROMOSOME MAINTENANCE 2) exhibited a high linkage disequilibrium (LD) with the respective SNPs. Further SNPs were found in an intergenic region of chromosome 2 between the two genes AT2G38630 and AT2G38640. All these SNPs had low minor allele frequencies (MAF) < 0.03, and thus represent rare SNPs in the studied set of accessions. AMM analysis of the O 3 /Munich full dataset (127 accessions) detected an additional SNP associated with the extent of O 3 -induced leaf damage in the 3' UTR of AT3G53400. However, GWAS with the KW test did not result in the identification of SNPs in either O 3 datasets (Supplementary Table S8).
GWAS of the full NO 2 dataset with 216 accessions did not reveal any signi fi cant SNPs using AMM analysis (Supplementary Table S8). However, when applying the KW test tens of SNPs in chromosome 1 showed significant association with the leaf damage score. The smallest P-value SNPs were localized in genes coding for a DNA glycosylase (AT1G19480), Transducin/WD40 repeat-like superfamily (AT1G19485), and a basic-leucine zipper transcription factor (AT1G19490). Other significantly phenotype-associated SNPs were in the coding sequences of AT2G03740 and AT4G32105 and on chromosome 5 between the genes AT5G08005 and AT5G08010.
In GWAS multiple testing correction is commonly used to provide fewer false positives (here Bonferroni adjustment). To allow comparisons between treatments, in subsequent analysis we also included SNPs with lower P-values. SNPs from the two O 3 experiments revealed no common SNPs with -log 10 P-values > 8, whereas 3 SNPs and 341 SNPs were shared with -log 10 P-values > 6 and > 4 ( Figure 5A, Supplementary Table S8). Altogether 47 (O 3 Helsinki) and 63 (O 3 Munich) SNPs showed -log 10 P-values > 6 and 1716 and 1681 SNPs showed -log 10 Pvalues > 4, respectively. Surprisingly, plotting the GWAS results of the NO 2 experiment against those of the O 3 experiments showed no shared SNPs with -log 10 P-values > 4 ( Figures 5B, C). In sum, the GWAS indicated very few genomic regions containing SNPs that have strong associations with O 3 or NO 2 leaf damage. Instead, there were many weak associations with O 3 and NO 2 sensitivity, and no overlap of SNPs between O 3 and NO 2 treatments.

Phenotyping of T-DNA insertion lines and O 3 sensitive mutants
We identified T-DNA insertion mutants for several GWAS candidate genes (for location of T-DNA inserts within the genes, see Supplementary Figure S3). We mainly selected based on Pvalues (Table 1), but for P-values with lower significance we also included candidate genes with biological functions related to plant defense responses. We measured the extent of damage as relative ion leakage after treatment for 6 h with 350 ppb O 3 or 10 ppm NO 2 (Table 1, Figure 6 and original data in Supplementary  Table S9). Relative ion leakage provides a better presentation of the data as the O 3 experiment was done in Helsinki and the NO 2 experiment in Munich. Col-0 was set as the baseline as all the mutants were in Col-0 background. Therefore, mutant line more sensitive than Col-0 have positive ion leakage values and more tolerant negative values.
Of the selected O 3 candidate genes, abh1 and cngc1 were O 3 sensitive (Supplementary Figure S4), whereas At3g53400 was NO 2 sensitive and At2g38640 more tolerant to NO 2 than the background accession Col-0. For mpk6 and five other mutants no significant change in ion leakage was observed compared to Col-0. Among the three NO 2 candidate genes tested, mutant lines with T-DNA insertions in AT1G19480 and AT1G19485 were O 3 and NO 2 sensitive, respectively. Next, we asked whether higher NO 2 concentrations would improve the detection of altered mutant phenotypes. Therefore, we exposed Col-0 to 30 ppm NO 2 which resulted in 50-65% ion leakage. With this treatment two tested NO 2 candidate were shown to be NO 2 sensitive (AT1G19480 and AT1G19485). While we could measure statistically significant changes in ion leakages for some of the mutants selected from GWAS candidates, overall these changes were small compared to the cell death found in O 3 or NO 2 sensitive mutants in Col-0 background (Brosche et al., 2014;Xu et al., 2015a;Kasten et al., 2016). This is in line with the overall GWAS results, i.e. O 3 and NO 2 damage in natural accessions are associated with many small effect loci rather than major effect loci.
Previously published O 3 sensitive mutants include ascorbate-deficient vtc1 (Conklin et al., 1996), ethylene overproducer eto1 (Rao et al., 2002), transcriptional coregulator rcd1 (Overmyer et al., 2000;Brosche et al., 2014), JA receptor coi1 (Xu et al., 2015a) and guard cell S-anion channel slac1 (Vahisalu et al., 2008). To allow a direct comparison of damage in mutants for GWAS candidates with the previously identified mutants, we used five O 3 sensitive mutants and performed ion leakage with 10 ppm NO 2 (we did not repeat the O 3 damage for these mutants, as they are extensively characterized in previous publications). The rcd1, vtc1 and slac1 mutants were strongly NO 2 sensitive, whereas eto1 and coi1 were not ( Figure 6). Kasten et al. (2016) tested 30 ppm NO 2 dose for these mutants and eto1 showed a NO 2 sensitivity phenotype whereas coi1 was similar to Col-0. In sum, the phenotyping of GWAS T-DNA mutant lines and O 3 sensitive mutants indicated that O 3 -and NO 2 -induced leaf damage was controlled by partially overlapping sets of genes. However, there were also exceptions, as the O 3 sensitive coi1 was tolerant to NO 2 .

Discussion
O 3 has emerged as a large threat to agricultural production in Asia, including wheat and rice (Feng et al., 2022b). Further understanding of the genetic basis of plant sensitivity to air pollutants is required to guide breeding programs aimed at providing plants with improved tolerance (Frei, 2015). O 3 tolerance and sensitivity traits are present in different genotypes and mapping populations of wheat (Feng et al., 2022a), maize (Choquette et al., 2019), rice (Frei et al., 2008) as well as A. thaliana [this study Jakobson et al., 2016;Morales et al., 2021)]. Combined genetic analysis with physiological traits has identified some of the mechanisms behind differential O 3 sensitivity; for example, a mutation that leads to more open stomata in the A. thaliana accession Cvi-0 leads to higher O 3 uptake and O 3 damage

FIGURE 5
Scatterplots for the comparison and correlations between -log 10 P-values of SNPs from genome-wide association analysis (AMM) between the three experiments with 119 common accessions. (A) plotted P-values from O 3 fumigations conducted in Munich (x-axis) and in Helsinki (y-axis), (B) plotted P-values from O 3 fumigations in Munich (x-axis) and NO 2 fumigations (y-axis), (C) plotted P-values from O 3 fumigations in Helsinki (x-axis) and NO 2 fumigations (y-axis). The areas with the smallest P-values (-log 10 P-value > 4) have coloured background. TABLE 1 Identified GWAS candidate genes, -log10 P-value of the most significant SNP in the regions, minor allele frequency and count of the most significant SNPs, significant differences in ion leakage measurements of the candidate gene mutants relative to Col-0 control after NO2 and O3 fumigations.  Jakobson et al., 2016). Variation in photosynthetic parameters offers another possibility to follow O 3 sensitivity traits (Choquette et al., 2019;Morales et al., 2021). Here we used the two air pollutants O 3 and NO 2 in A. thaliana to further understand their impact on regulation of signaling pathways and PCD. Ultimately, this could provide new mechanisms for tolerance that could be targeted in plant breeding programs.

O 3 -versus NO 2 -induced transcriptional responses
Exposure of Col-0 to 10 ppm NO 2 for 1 h or 350 ppb O 3 for 2 h resulted in massive transcriptional changes (Xu et al., 2015b;Mayer et al., 2018). Both gases regulated largely overlapping sets of genes involved in ROS, ethylene signaling, pathogen resistance, and cell death ( Figure 1B, Supplementary Tables S5  and S6). This is consistent with the rapid accumulation of ROS, NO, and ethylene in O 3 -exposed tobacco (Ederli et al., 2006). ROS and RNS bursts were also observed after NO 2 fumigation of A. thaliana (Kasten et al., 2016). Plants produce ROS and RNS molecules as signaling molecules to regulate local and systemic long distance defence responses (Wang et al., 2013;Waszczak et al., 2018;Hancock and Neill, 2019;He et al., 2022). Accordingly, O 3 and NO 2 both increase transcript levels for pathogen responsive genes (Xu et al., 2015a;Mayer et al., 2018). In sum, O 3 and NO 2 probably act both as donors (i.e. to generate) as well as inducers of simultaneous ROS and RNS bursts that ultimately lead to the onset of defence responses.
We focused on the GO category cell death and performed Bayesian hierarchical clustering with O 3 and various NO related treatments ( Figure 1A). This revealed similar expression profiles in both O 3 tolerant (Col-0, C24) and sensitive accessions (Cvi-0, Te). To corroborate this finding, we analysed transcript levels for five cell death related marker genes in a side-by-side comparison of nine accessions treated with 350 ppb O 3 or 10 ppm NO 2 for 2 h (Figure 2). Both gases increased transcript levels of FMO1, GRX480, CEJ1, and RAP2.6 that function in defence signaling (Krishnaswamy et al., 2011;Caarls et al., 2015;Hartmann et al., 2018). Importantly, RBOHF was differentially regulated with increased transcript levels after O 3 , but decreased after NO 2 treatment ( Figure 2). RBOHF, but not RBOHD, was previously implicated as a regulator of O 3 cell death (Xu et al., 2015a), hence differential use of ROS produced from RBOHs could be a mechanism to regulate cell death in response to different signals. The transcript level variation between the accessions was independent of their O 3 or NO 2 sensitivity. For instance, the O 3 tolerant accession Ts-1 showed comparable gene regulation Percentage of change in ion leakage in GWAS candidate mutants (the upper part) and known O 3 sensitive mutants (the lower part of the figure) relative to Col-0. Treatment with 350 ppb of O 3 for 6h is marked in turquoise, 10 ppm of NO 2 for 6h in orange and 30 ppm of NO 2 for 1h in gray. Statistical significance is marked with asterisks: ***P < 0.001, **P < 0.01 and *P < 0.05.
as the sensitive accession Cvi-0. Taken together, Figures 1, 2 support the conclusion that O 3 and NO 2 regulate the expression of largely overlapping sets of defence-related genes also in different A. thaliana accessions. This implies that the mechanism(s) used by A. thaliana to perceive ROS (O 3 ) and RNS (NO 2 ) are conserved in genetically distant A. thaliana accessions. Initial NO perception in A. thaliana takes place via targeted degradation of group VII ethylene response factors (ERFs) (Gibbs et al., 2014). Further down-stream signaling is proposed to be mediated by several other transcription factors, including RAP2.6 (Imran et al., 2018;Leon et al., 2020). Since RAP2.6 transcript levels were enhanced by both O 3 and NO 2 , it represents a target for both ROS and RNS signaling.
There was one exception to the common transcriptional response by O 3 and NO 2 , RBOHF was differentially regulated, with increased transcript levels after O 3 , but decreased after NO 2 treatment ( Figure 2). In both plant stress, PCD and developmental responses, the RBOH proteins produce superoxide as signaling molecules (Castro et al., 2021). Out of the ten A. thaliana RBOHs, RBOHD and RBOHF are the main producers of ROS during various aspects of defence signaling (Castro et al., 2021). RBOHF was previously implicated as a regulator of O 3 cell death (Xu et al., 2015a), and in ROS transcriptional responses (Willems et al., 2016). ROS from RBOHD and RBOHF are also required to execute cell death in pathogen HR and progression of cell death lesions (Torres et al., 2002;Torres et al., 2005). Differential use of ROS produced from distinct RBOHs could be a mechanism to regulate cell death in response to different signals. The opposite regulation RBOHF transcript levels were present in all tested accessions (Figure 2), this means that A. thaliana can activate distinct signaling pathways from O 3 versus NO 2 . Further analysis of the promoter region of RBOHF could lead to identification of e.g. transcription factor(s) and promoter elements that regulate this specific signaling pathway with contrasting regulation by O 3 and NO 2 .

Use of GWAS to identify genes regulating O 3 and NO 2 cell death
Improvement of plant varieties by breeding could help to reduce yield losses due to pollutant-induced leaf damage. Plant breeding can be guided by the results from quantitative trait loci (QTL) or association mapping studies (Frei, 2015;Ueda et al., 2015;Begum et al., 2020). Here, we used A. thaliana natural accessions to identify genes involved in regulation of lesion formation after O 3 and NO 2 exposure. We performed GWAS on two independent O 3 screens and one NO 2 screen. 119 accessions were investigated in all three experiments. Initial GWAS runs focused only on 119 accessions assuming that at least some genetic loci would be associated with both O 3 -as well as NO 2induced leaf phenotypes. However, this approach did not reveal any significant SNPs shared between the O 3 and NO 2 screens (Figure 4).
Even with standardized growth conditions, growth and molecular responses of A. thaliana show differences between different laboratories (Massonnet et al., 2010). To explore these differences, we performed the phenotyping of O 3 tolerance at two different facilities (Helsinki and Munich). Replication of GWAS is common in human studies but relatively rare in other organisms. As we repeated our O 3 GWAS at two different facilities, this can give some insight into what to expect from GWAS replications in A. thaliana. There was some phenotypic variation, but also a significant positive correlation between the same accessions in the two O 3 experiments. Importantly, we found noticeable overlap in small P-value GWAS SNPs between the two O 3 datasets ( Figure 5A, 341 shared SNPs with -log 10 Pvalues > 4), which emphasizes that A. thaliana has a robust genetic response to O 3 . For instance, the most significant SNP for the common 119 accessions in Helsinki (AT3G61410) had also a small P-value in Munich (-log 10 P-value = 4.1). Furthermore, by combining results from shared small P-value SNPs from GWAS replications, additional candidate genes can be selected for further study. This would be especially useful as many traits of interest are controlled by several small effect genes, for which SNPs may not reach significance level after correction for multiple testing.
Subsequently, we analysed the datasets individually and we found 12 genomic loci significantly associated with the extent of leaf damage ( Figure 5 and Table 1) but again none of these loci was shared between the NO 2 and O 3 screens. GWAS for O 3 tolerance has previously been performed in rice (Ueda et al., 2015), which identified 16 loci with rather weak phenotypic associations (P<0.0001). Similarly, in a study with 150 wheat varieties statistics indicated weak associations between SNPs and O 3 -induced phenotypes because the determined P-values were rarely below 0.0001 (Begum et al., 2020). Classical QTL mapping studies indicated that several genes in different chromosomal locations control O 3 sensitivity in A. thaliana Xu et al., 2015b;Jakobson et al., 2016). These findings argue for O 3 and NO 2 phenotypes being determined by multiple small-effect loci that are detectable either by QTL-mapping families with contrasting parental phenotypes (which gives the possibility to identify rare alleles in A. thaliana populations) or by using large numbers of accessions in GWAS to improve statistical sensitivity (and more likely to identify common alleles). Accordingly, GWAS with O 3 screening data from 127 accessions (Munich screen) resulted in the identification of only a single significantly associated genomic region (only SNPs in one gene in high linkage disequilibrium with the most significant SNP) whereas data from 372 accessions (Helsinki screen) revealed 3 associated genomic regions (with SNPs in 8 genes in high LD with the most significant SNPs). Hence, future screens should include as many of the >1000 sequenced accessions (Alonso-Blanco et al., 2016), as possible to identify more genes linked to O 3 -and NO 2 -associated phenotypes and assess whether there is reasonable overlap in the genetics underlying responses to both pollutants. However, with the current datasets we could not find significantly, or even among small P-value ( Figures 5B, C), overlapping SNPs and thus it appears that O 3 versus NO 2 induced lesions are controlled by different genetic loci.
Based on the GWAS results we screened 13 T-DNA knockout mutants for candidate genes (Table 1 and Figure 6). Five mutants were sensitive and one mutant tolerant as compared to wild-type plants. ABH1 is a gene coding for a nuclear mRNA cap-binding protein that participates in abscisic acid signaling and mRNA processing (Hugouvieux et al., 2001;Laubinger et al., 2008). CNGC1 is an ion channel likely functioning in calcium signaling (Sunkar et al., 2000). The mutant for At1g19480 was sensitive to both 10 ppb O 3 as well as 30 ppm NO 2 ( Figure 6). AT1G19480 might be involved in DNA repair that is an important process in ROS-exposed plants (Nisa et al., 2019). CPuORF46 (AT3G53400) encode an upstream open reading frame (uORF) which can conditionally regulate translation of the main ORF, and has been shown to be responsive to heat stress (Causier et al., 2022). Future research can characterize the physiological mechanisms causing the NO 2 -and O 3 -phenotypes in these mutants including measurements of the stomatal aperture, transcriptional changes, stress hormone levels, and antioxidant status.
Both O 3 and NO 2 cause formation of PCD lesions, with similarities to pathogen HR (Overmyer et al., 2005;Kasten et al., 2016). To further compare how O 3 and NO 2 regulate lesion formation we used the O 3 sensitive vtc1 and slac1 (Conklin et al., 1996;Vahisalu et al., 2008). Both mutants displayed stronger cell death and increased ion leakage upon exposure to NO 2 than any other mutant tested in the current study ( Figure 6B). The vtc1 mutant is ascorbate-deficient whereas slac1 has an increased stomatal aperture and impaired responses to signals leading to stomatal closure. Hence, antioxidant levels and stomatal regulation restricting entry of air pollutants into the plant are important common determinants of both O 3 and NO 2 toxicity.
Overall, using O 3 and NO 2 we show that ROS and RNS have largely overlapping transcriptional responses, but at the same time, they also have distinct signaling roles as exemplified by the contrasting transcriptional regulation of RBOHF. Similarly, mutant analysis also revealed mutants that were sensitive to both gasses as well as mutants sensitive to only one gas. No common SNPs were identified by the GWAS analysis suggesting that natural variation in the strength of PCD induced by the Signalling events initiated by O 3 and NO 2 . O 3 and NO 2 enter the plant through stomata, and in the apoplast break down into ROS (including H 2 O 2 and O 2 .-), and RNS (including NO). This initiates various signalling pathways leading to large-scale transcriptional changes, where O 3 and NO 2 mostly regulate the same set of genes ( Figure 1). In sensitive genotypes, PCD is activated leading to leaf damage. Whether PCD is activated depends on the balance of ROS and RNS signalling (Delledonne et al., 2001;Ahlfors et al., 2009), the antioxidant capacity to remove ROS and RNS (Conklin et al., 1996;Kasten et al., 2016), hormone signalling (Xu et al., 2015a;Kasten et al., 2016) and other to be discovered mechanisms. Placement of new candidate PCD regulators identified from GWAS into signalling pathways will require additional studies. The application of O 3 and NO 2 leads to a clearly defined subcellular source of ROS and RNS in the apoplast. Multiple other abiotic and biotic stresses also lead to active ROS and RNS signalling (Wendehenne et al., 2014;Waszczak et al., 2018;He et al., 2022;Hussain et al., 2022). However, in response to stress the production site of ROS and RNS can come from multiple parts of the cell including apoplast, chloroplast and mitochondria (Waszczak et al., 2018;Hancock and Neill, 2019;Castro et al., 2021). Created with BioRender.com. gasses is caused by different small effect genes. Due to the relative simplicity of applying O 3 or NO 2 to plants in controlled growth conditions, further studies with these gasses will allow dissection of ROS and RNS signaling pathways in plant stress and PCD regulation. Identification of novel PCD regulators can aid in development of strategies to combat the negative effects of air pollution. Figure 7 illustrates connections between O 3 and NO 2 and other stresses that can initiate defence signaling and PCD.

Data availability statement
The original contributions presented in the study are included in the article and supplementary material. Raw data used for re-analysis of microarray and RNA-seq data are described in detail in Materials and Methods and in Supplementary Table S2.

Funding
For funding, in Helsinki, we acknowledge a grant number 307335 from the Academy of Finland, for Centre of Excellence in Molecular Biology of Primary Producers. The gas exposure screens in Munich were supported by the European Plant Phenotyping Network (EPPN, Project No. 28443).