Distinct transcriptional profiles of ozone stress in soybean (Glycine max) flowers and pods

Background Tropospheric ozone (O3) is a secondary air pollutant and anthropogenic greenhouse gas. Concentrations of tropospheric O3 ([O3] have more than doubled since the Industrial Revolution, and are high enough to damage plant productivity. Soybean (Glycine max L. Merr.) is the world’s most important legume crop and is sensitive to O3. Current ground-level [O3] are estimated to reduce global soybean yields by 6% to 16%. In order to understand transcriptional mechanisms of yield loss in soybean, we examined the transcriptome of soybean flower and pod tissues exposed to elevated [O3] using RNA-Sequencing. Results Elevated [O3] elicited a strong transcriptional response in flower and pod tissues, with increased expression of genes involved in signaling in both tissues. Flower tissues also responded to elevated [O3] by increasing expression of genes encoding matrix metalloproteinases (MMPs). MMPs are zinc- and calcium-dependent endopeptidases that have roles in programmed cell death, senescence and stress response in plants. Pod tissues responded to elevated [O3] by increasing expression of xyloglucan endotransglucosylase/hydrolase genes, which may be involved with increased pod dehiscence in elevated [O3]. Conclusions This study established that gene expression in reproductive tissues of soybean are impacted by elevated [O3], and flowers and pods have distinct transcriptomic responses to elevated [O3]. Electronic supplementary material The online version of this article (doi:10.1186/s12870-014-0335-y) contains supplementary material, which is available to authorized users.


Background
Current tropospheric O 3 concentrations ([O 3 ]) are estimated to cost $14 to $26 billion in annual global crop economic losses [1] and severely impact human health, accounting for an estimated 0.7 million deaths per year [2]. Ozone in the troposphere is formed through the photochemical oxidation of volatile organic compounds (VOCs), carbon monoxide and methane in the presence of nitrogen oxides (NO x ) [3]. Ozone is a dynamic pollutant and concentrations vary temporally and spatially, with higher concentrations in the Northern Hemisphere compared to the Southern Hemisphere, and typically higher [O 3 ] in the summer compared to the winter [3]. Background tropospheric [O 3 ] have more than doubled since the Industrial Revolution and are projected to increase by an additional~20% by the year 2100 if current high emission rates continue [4]. In the crop growing regions of the Northern Hemisphere, summer concentrations of O 3 often exceed 40 ppb, which exceeds the critical threshold for damage to sensitive crops, including soybean (Glycine max) [5].
When taken up by plants, O 3 is converted into other reactive oxygen species (ROS), and can induce signaling pathways that lead to programmed cell death, especially with exposure to very high [O 3 ] [6]. At lower concentrations, chronic exposure to elevated [O 3 ] decreases photosynthetic carbon assimilation and stomatal conductance, and accelerates the process of senescence [7,8]. In addition to leaf-level effects, O 3 negatively impacts plant fitness and reproductive development, which can be mediated through reduced carbon allocation from source tissues and/or through direct effects on reproductive tissues [9,10]. A meta-analysis of published studies from 1968 to 2010 of O 3 effects on plant reproductive processes reported that exposure to elevated [O 3 ] decreased seed number and seed size, as well as fruit number and fruit size when compared to plants grown in charcoal-filtered, O 3 -free air [11]. However, the meta-analysis also showed that elevated [O 3 ] did not significantly alter inflorescence number, flower weight or flower number [11]. This suggests that plants can compensate to some extent from O 3 damage [12], and also that the effects of O 3 can be tissue-specific.
Soybeans have naturally high levels of floral and pod loss, and subsequent seed and yield loss is greatest when stress occurs during flower and early pod development [13]. Flower and pod abscission can range from 32 to 82% in soybean [14][15][16], but this varies considerably with location on the plant [15][16][17], location in the canopy [18], source-sink relations [19], hormone levels [13,20,21], shade [22] and water status [13,23,24]. Ethylene promotes flower and pod abscission in soybean [25], and elevated [O 3 ] can increase ethylene emission in plants [26]. Therefore, elevated [O 3 ] has the potential to increase flower and pod abscission. In field-grown soybean exposed to elevated [O 3 ] for an entire growing season [27], pod production was decreased by elevated [O 3 ], but flower number was not affected (Figure 1). Based on this evidence from the field, it is hypothesized that the

Results and discussion
Overlapping effects of elevated [O 3 ] on the transcriptome of flower and pod tissue in soybean Flower and pod development in soybean are sensitive to environmental stress [23,24,38,39], and elevated [O 3 ] significantly impacted pod production, but not flower production ( Figure 1). In order to identify the genetic mechanisms underpinning O 3 response in soybean pods and flowers, the transcriptome of flower and pod tissues was compared using RNA-Sequencing (RNA-Seq). The global transcriptional analysis showed the magnitude of potential responses to elevated [O 3 ] in flowers and pods was similar, with genes showing approximately the same range of both mean expression values in flowers and pods, and similar potential log fold change responses to elevated [O 3 ] in the two tissues ( Figure 2). However, more than three times as many genes were differentially expressed in flower tissue (4,595 genes) than in pod tissue (1,375 genes; Figure 3    Many of the transcripts that fell on the 1:1 line in Figure 5 were involved in signaling and RNA processing, including 12 leucine-rich repeat receptor-like kinases (RLKs) and 3 cysteine-rich Domain of Unknown Function 26 (DUF26) RLKs (also known as cysteine-rich receptor-like kinases, CRK). Plant RLKs are transmembrane proteins involved in signal perception and form a large multi-gene family with regulatory roles in development, abiotic and biotic stress responses in plants [40,41]. Recent analysis of the response of Arabidopsis DUF26 RLKs showed that many of the 44 RLKs were specifically up-regulated in response to O 3 stress in leaves [42], including DUF26 30 (CRK 26), DUF26 29 (CRK 29) and DUF26 41 (CRK 2), which also had a significant increase in expression in soybean pods and flowers exposed to elevated [O 3 ]. Wraczek et al. [42] found that the general pattern of DUF26 expression responses to O 3 was most similar to the transcriptional response to pathogen infection, which like O 3 elicits an ROS burst in the apolost. The transcriptional response to O 3 however, was very different from expression responses to high light treatments or chemical treatments that increased ROS production in chloroplasts or mitochondria [42]. Thus, it was further suggested that the DUF26 domain, which has a conserved cysteine motif C-8X-C-2X-C, could act as an apolastic ROS sensor [42]. A number of WRKY domain transcription factors also showed significantly greater expression under elevated [O 3 ] in both soybean flowers and pods (Glyma04g40130, Glyma06g14720, Glyma14g36430, Glyma14g36438, Gly-ma14g36446). The WRKY transcription factor family is one of the largest families of transcription factors in plants, with 133 members in the soybean genome [43]. WRKYs function in many plant processes including response to biotic and abiotic stresses, and senescence [44]. Up-regulation of WRKY transcription factors in response to O 3 stress has been previous reported, primarily in the leaves of trees [45][46][47]. Two of the WRKY transcription factors with increased expression in both flower and pod tissues in response to elevated [O 3 ] (Gly-ma04g40130 and Glyma06g14720) were likely formed through a segmental duplication event that is estimated to have occurred 20 million years ago [43].

Distinct effects of elevated [O 3 ] on the transcriptome of flower and pod tissue in soybean
Although there was some overlap in transcriptional responses to elevated [O 3 ] in flowers and pods, the vast majority of genes changing in either tissue were distinct (Figure 3), and even among the genes that were expressed in both flower and pod tissues, the fold changes in expression were of different magnitudes or in the opposite direction ( Figure 5). In flower tissue, the genes with the greatest increase in abundance in response to elevated [O 3 ] included matrix metalloproteinase (MMP) genes and genes related to hormone metabolism and signaling (Table 1). Genes annotated as MMPs also had high mean expression levels (Figure 2). While little is known about their role in soybean flowers, in other tissues MMPs function in degradation of the extracellular matrix (ECM) in response to senescence, stress and programmed cell death [48][49][50][51]. Domain analysis indicated that 7 of the 9 differentially expressed genes annotated as MMPs had both a cysteine switch domain and a zinc-binding domain, both of which are required for characterization as a MMP ( Figure 6) [52][53][54]. Those genes with both required domains were termed putative soybean flower MMP genes.
The putative MMP gene Glyma02g03301 had two identical cysteine switch domains and zinc-binding domains, which was unique compared to the other putative flower MMP genes. All putative flower MMP genes had a signal peptide and transmembrane domain, with the exception of Glyma02g03250 (Figure 6; Additional file 1). Several of the putative flower MMP genes contained a GPI-anchor modification site, which was similar to 3 Arabidopsis MMP genes, At2-MMP, At4-MMP, At5-MMP [47], and the known soybean MMP gene GmMMP2 (referred to as Gm2-MMP) [42] (Figure 6). None of the putative flower MMP genes or known soybean genes (Gm2-MMP or SMEP1) [48,55,56] contained a furin cleavage site, which was present in several Arabidopsis MMP genes. When amino acid sequence similarity identity was compared between all putative flower MMP genes and Gm2-MMP using Clustal W (http://npsa-pbil.ibcp.fr/cgi-bin/npsa_automat.pl? page=/NPSA/npsa_server.html), little homology between the flower MMP and leaf MMP genes was found, with the exception of Glyma01g04350 which showed 99% sequence similarity to Gm2-MMP (data not shown). All putative flower MMP genes, with the exception of Glyma01g04350, had an E (glutamate) to Q (glutamine) residue substitution in the zinc-binding motif of the catalytic domain, which has been identified in other legume species [51]. The glutamate residue is required for functional protease activity [57], thus the amino acid switch in soybean flower MMPs may render these inactive. Still, they may be important for O 3 stress response because experiments with Medicago truncatula have also demonstrated a functional role for proteolytically-inactive MMPs in biotic stress response [51].
The responsiveness of putative soybean flower MMPs to elevated [O 3 ] is consistent with the ECM being the primary point of O 3 contact within plant cells and the location where antioxidant metabolism begins to protect cells from ROS damage [34,58]. Stress-responsive signaling pathways, including jasmonic acid, salicylic acid and ethylene-dependent redox signaling are all triggered by the redox sensing that occurs in the ECM [34]. The soybean MMP gene Gm2-MMP was up-regulated consistently with the release of ROS during pathogenic infection [50], possibly linking ROS signaling and MMP gene expression in soybean stress response. Previous analysis of Arabidopsis MMP gene expression revealed that At3-MMP was expressed at greater abundance in response to O 3 treatment, with a slight increase in At2-MMP in response to O 3 as well [59]. While the putative MMP genes are present in high abundance in flower tissues exposed to elevated [O 3 ], analysis of the expression profiles of the putative MMP genes in soybean using RNA-Seq Atlas (http://soybase.org/soyseq/) found that these genes were not present, or present in low abundance in other soybean tissues. Therefore, it is hypothesized that the increase in abundance of the putative MMP genes identified in this study may represent a distinct flower response to O 3 stress in soybean.
In pod tissue, cell wall modification and calcium signaling genes showed the greatest increase in abundance in response to elevated [O 3 ] ( Table 2). Gene ontology (GO) enrichment analysis of biological processes was performed for genes differentially expressed only in pod tissue. Apoptosis, signal transduction, ATP biosynthetic processes, cellular glucan metabolic processes, protein amino acid phosphorylation and innate immune responses were enriched in pod tissue (Additional file 2). These activities are known to increase in plants in response to both abiotic and biotic stress [60][61][62][63][64][65], and the possibility that O 3 stress co-opts pathways involved in biotic stress response has been previously proposed [66,67]. The genes with the greatest increase in abundance in response to elevated [O 3 ] were xyloglucan endotransglucosylase/hydrolase family proteins (XTH) ( Table 2). Genes annotated as XTHs also had high mean expression, along with the greatest increase in abundance in response to O 3 in pod tissue (Figure 2). These genes belong to the GO biological process of cellular glucan metabolic processes, which is highly enriched in pod tissues (Figure 7). Analysis of the putative XTH genes in soybean using RNA-Seq Atlas (http://soybase. org/soyseq/) showed that these genes were not present or in low abundance in other tissues, indicating that these genes may represent a distinct pod response to elevated [O 3 ]. XTH is a well-known cell wall-modifying enzyme that plays a role in growth and differentiation in plants [68]. Genes in the XTH family are involved in cell elongation in vascular cells [69,70], epidermal cells [71], inflorescence apices [72], primary roots [73] and during somatic embryogenesis [74]. XTH also plays a role in floral organ abscission [75,76]. Due to the similarity of flower abscission and pod dehiscence zone [77] and the known response of XTH genes to oxidative [78], water [79] and biotic stress [80], it is hypothesized that XTH genes may play a role in pod dehiscence in soybean exposed to elevated [O 3 ].

Conclusion
Soybean is an O 3 -sensitive crop, with current tropospheric [O 3 ] costing billions of dollars in lost production annually. In this study, it was established that gene expression in reproductive tissues in soybean is altered by elevated [O 3 ]. There were 4,703 transcripts responsive to elevated [O 3 ] in both flower and pod tissues, yet those genes did not respond consistently in the two tissues. This indicates that reproductive tissues have more distinct than similar transcriptomic responses to elevated [O 3 ]. Flower tissues responded to elevated [O 3 ] through increased expression of MMP genes. It was notable that these flower MMP genes may not be proteolytically active based on amino acid composition, but they clearly respond to O 3 stress. Pod tissues responded to elevated [O 3 ] through increased expression of cell expansion genes. The increased transcript abundance of XTH genes supports a role of these genes in pod dehiscence in soybean exposed to elevated [O 3 ].

Growth chamber experimental design and conditions
Soybean (Glycine max L. Merr. cv. 93B15; Pioneer Hi-Breed) was grown in ambient air (<20 ppb) and elevated ozone (150 ppb) in 14 h/10 h day/night schedules under PPFD of~650-750 μmol m −2 s −1 ; RH 60%; 25°C day/21°C night conditions in 8 growth chambers (Conviron, Winnipeg, Manitoba, Canada). Soybean plants were grown in 6-L pots (Classic C600, Nursery Supplies, Chambersburg PA, USA) in sterilized soil (LC-1 Sunshine Mix (SunGro Horticulture Canada Ltd, Bellevue, WA, USA)) and treated with 50% Long Ashton solution supplemented with 3 mM NH 4 NO 3 [81]. Two seeds were planted per pot~4 cm below the soil surface and then thinned to one plant per pot once seeds successfully germinated. A total of 12 plants were grown per chamber in a randomized complete block design (n = 4).
Plants were rotated among chambers once a week and within chambers every two days to minimize chamber effects.

Tissue sampling and molecular analyses
Tissue sampling for RNA was done during R2 (full bloom) and R4 (full pod) for growth chamber grown plants. Plants were considered at full bloom when there was an open flower at one of the first two uppermost nodes with a fully expanded leaf. Plants were considered at full pod when there was a pod 2 cm in length present on one of the four uppermost nodes with a fully expanded leaf. At each stage the appropriate tissue was sampled (full open flowers at R2 and initiating pods at R4). Sampling was done at the nodes 2-4 (from the top of the plant) in order to avoid compensation and senescence effects on the upper and lowermost nodes. Tissue from four plants was sampled per developmental stage per block. Immediately after collection flower and pod tissue was plunged into liquid N and stored at −80°C. Flower or pod tissue was ground to a fine power using a mortar and pestle. Total RNA was extracted from ground tissue using Pure-Link Plant RNA Reagent (Ambion, by Life Technologies Corp., Grand Island, NY, USA) according to the manufacturer's protocol. RNA quantity was determined with a spectrophotometer (Nanodrop 1000, Thermo Fischer Scientific, Waltham, MA, USA) and RNA quality was Analysis of gene ontology (GO) term enrichment of biological processes containing XTH genes in pod tissues. Biological terms with increasing overrepresentation in pod tissues exposed to elevated [O 3 ] are represented by increasingly red colors. GO term enrichment was performed using single enrichment analysis (SEA) tool on AgriGo (http://bioinfo.cau.edu.cn/agriGO/). assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) on an RNA Nano chip. Genomic DNA contamination was removed from RNA samples using Turbo DNase treatment (Applied Biosystems/Ambion, Austin, TX, USA) according to the manufacturer's protocol. cDNA library preparation was done using the Illumnina TruSeq Sample Prep kits (Illumina Inc. San Diego, CA, USA). Each library fragment was barcoded during library preparation and multiplexed for sequencing. Tissue samples per block (4 subsamples) were pooled for a total of 8 libraries prepared for each tissue (16 libraries total).

RNA-sequencing (RNA-seq), bioinformatics and statistical analysis
Sequencing was done at the Roy J. Carver Biotechnology Center using the Illumina Genome HiSeq 2000 (Illumina Inc. San Diego, CA, USA, http://www.illumina.com) and Cassava pipeline 1.8 to obtain 100 nt single-end reads. Samples were sequenced in groups of 4 across 4 lanes and generated~31-63 million reads per sample. All FASTQ files from all sequencing runs are located in the Small Read Archive (http://www.ncbi.nlm.nih.gov/sra), SRP035871, BioProject number PRJNA236472. Quality control for reads generated from sequencing was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/). Sequenced reads were aligned to the soybean reference genome (Gmax_189.fa, www.phytozome. net) using Bowtie [82]. All valid alignments per read were reported allowing up to three mismatches. Alignment summary statistics are presented in Additional file 3. Aligned sequence reads and a list of genomic features (Gmax_189_gene.gff3, www.phytozome.net) were input into HTSeq to generate read counts using the htseq-count and -m union options. These counts were then input into SAS (SAS Institute, Version 9.2, Cary, NC, USA) for normalization and statistical analysis. Genes with counts of 10 or less were removed from all subsequent statistical analyses. Read counts were normalized using the natural log of the upper quartile (ln_uq) [83,84]. All count data can be found in Additional files 4 and 5. Differential gene expression was determined using a mixed effects linear model Y ijkl = m + t i + γ j + ρ k + ɛ ijkl . Y is the normalized estimate of the expression for the fixed effect of condition (i = ozone/ambient), the random effect of block (j = 1,2,3,4) and the random effect of lane (k = 1,2,3,4). A log fold change represents the difference of the ln_uq normalized count data for elevated [O 3 ] minus ambient [O 3 ]. The assumptions of normality were tested using the Shapiro-Wilk test [85] for each gene. A multiple test correction was applied using the linear step-up method of [86]. Analyses were conducted in SAS (SAS Institute, Version 9.2, Cary, NC, USA). the cDNA library preparation. EAA designed the experiments, provided critical input in the transcriptomic analysis and wrote the manuscript. All authors read and approved the final manuscript.