# Genome-Wide Transcriptional Profiling and Metabolic Analysis Uncover Multiple Molecular Responses of the Grass Species Lolium perenne Under Low-Intensity Xenobiotic Stress

Anne-Antonella Serra, Ivan Couée, David Heijnen, Sophie Michon-Coudouel, Cécile Sulmon, Gwenola Gouesbet

<figure><img src="image_1.png" /></figure>

To cite this version:

Anne-Antonella Serra, Ivan Couée, David Heijnen, Sophie Michon-Coudouel, Cécile Sulmon, et al. Genome-Wide Transcriptional Profiling and Metabolic Analysis Uncover Multiple Molecular Responses of the Grass Species Lolium perenne Under Low-Intensity Xenobiotic Stress. Frontiers in Plant Science, 2015, 6, pp.1124. 10.3389/fpls.2015.01124. hal-01245729

HAL Id: hal-01245729

https://univ-rennes.hal.science/hal-01245729

Submitted on 10 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Genome-Wide Transcriptional Profiling and Metabolic Analysis Uncover Multiple Molecular Responses of the Grass Species *Lolium perenne* Under Low-Intensity Xenobiotic Stress

Anne-Antonella Serra<fnref n="1" />, Ivan Couée<fnref n="1" />, David Heijnen<fnref n="1" />, Sophie Michon-Coudouel², Cécile Sulmon<fnref n="1" /> and Gwenola Gouesbet<fnref n="1" />*

## OPEN ACCESS

**Edited by:**
Girdhar Kumar Pandey,
Delhi University South Campus, India

**Reviewed by:**
Eric Van Der Graaff,
Copenhagen University, Denmark
Manoj Prasad,
National Institute of Plant Genome
Research, India

*Correspondence:
Gwenola Gouesbet
gwenola.gouesbet@univ-rennes1.fr

**Specialty section:**
This article was submitted to
Plant Physiology,
a section of the journal
Frontiers in Plant Science

**Received:** 02 June 2015
**Accepted:** 27 November 2015
**Published:** 17 December 2015

**Citation:**
Serra A-A, Couée I, Heijnen D, Michon-Coudouel S, Sulmon C and Gouesbet G (2015) Genome-Wide Transcriptional Profiling and Metabolic Analysis Uncover Multiple Molecular Responses of the Grass Species *Lolium perenne* Under Low-Intensity Xenobiotic Stress.
Front. Plant Sci. 6: 1124.
doi: 10.3389/fpls.2015.01124

*Lolium perenne*, which is a major component of pastures, lawns, and grass strips, can be exposed to xenobiotic stresses due to diffuse and residual contaminations of soil. *L. perenne* was recently shown to undergo metabolic adjustments in response to sub-toxic levels of xenobiotics. To gain insight in such chemical stress responses, a *de novo* transcriptome analysis was carried out on leaves from plants subjected at the root level to low levels of xenobiotics, glyphosate, tebuconazole, and a combination of the two, leading to no adverse physiological effect. Chemical treatments influenced significantly the relative proportions of functional categories and of transcripts related to carbohydrate processes, to signaling, to protein-kinase cascades, such as Serine/Threonine-protein kinases, to transcriptional regulations, to responses to abiotic or biotic stimuli and to responses to phytohormones. Transcriptomics-based expressions of genes encoding different types of SNF1 (sucrose non-fermenting 1)-related kinases involved in sugar and stress signaling or encoding key metabolic enzymes were in line with specific qRT-PCR analysis or with the important metabolic and regulatory changes revealed by metabolomic analysis. The effects of pesticide treatments on metabolites and gene expression strongly suggest that pesticides at low levels, as single molecule or as mixture, affect cell signaling and functioning even in the absence of major physiological impact. This global analysis of *L. perenne* therefore highlighted the interactions between molecular regulation of responses to xenobiotics, and also carbohydrate dynamics, energy dysfunction, phytohormones and calcium signaling.

**Keywords:** RNA-Seq, glyphosate, tebuconazole, no observable adverse effect, residual pollution, perennial ryegrass, SnRKs

**Abbreviations:** G, glyphosate; T, tebuconazole; GT, combination of glyphosate and tebuconazole.

<footnote n="1">¹Centre National de la Recherche Scientifique, Université de Rennes 1, UMR 6553 ECOBIO, Rennes, France, ²Centre National de la Recherche Scientifique, Université de Rennes 1, UMS 3343 OSUR, Rennes, France</footnote>
# INTRODUCTION

Modern agriculture uses large amounts of phytosanitary products to maximize crop production. These products are mainly pesticides or mixtures of pesticides, which have been developed to repel, attenuate or kill pests and competitive plants (Arias-Estévez et al., 2008). However, runoff, leaching, or spray drift lead to the entry of a large fraction of pesticides into environmental compartments (soil, water, sediment, or atmosphere). These diffuse and residual contaminations are composed of mixtures of parent compounds, of their degradation products, and also of associated adjuvants (Helander et al., 2012) and induce many environmental impacts (Patty et al., 1997; Köhler and Triebskorn, 2013). Pesticides and related-degradation products chemically stress many non-target organisms in natural ecosystems, among which plants are particularly affected as sessile organisms (Child et al., 1993; Serra et al., 2013, 2015).

Mechanisms and regulations of chemical stress responses may differ between species as a result of micro- and macroevolutionary processes (Medina et al., 2007) and may influence the sensitivity or tolerance of plants and thereby their capabilities to maintain growth and development in polluted areas. Moreover, herbicide efficiency is increasingly affected by the emergence of resistance processes (Délye, 2013). Target-site resistance (TSR) results from gene mutations that alter herbicide targets, whereas processes of non-target-site resistance (NTSR) can result from multiple mechanisms of detoxification, tolerance or regulation (Délye, 2013). In the last decades, large-scale pesticide use on agricultural lands, resulting in long-term contamination (Singh et al., 2004), has put strong ecological and evolutionary pressure on the dynamics of plant communities, leading to resistance emergence (Délye, 2013). Widespread agricultural weeds are controlled by herbicide applications, which result in recurrent exposition for the surrounding vegetation generally composed of a majority of grasses. Among these species, several populations of *Lolium perenne* have been described as resistant to herbicides. Recently, it has been demonstrated that, in a population that displays glyphosate resistance, other mechanisms than mutation in the target site of glyphosate, the plasticid enzyme 5-enolpyruvylshikimate-3-phosphate synthase (EPSPS), were involved (Salas et al., 2015). Similarly, resistance to the acetolactate-synthase (ALS) inhibiting herbicide pyroxsulam in *Lolium* sp. populations involves a NTSR response implying differential gene expression and different mechanisms that remain to be elucidated (Duhoux et al., 2015). Multiple-herbicide resistance has also been described in a specific population of *Lolium perenne* spp. *multiflorum*, but although some resistances to chemical groups of ALS inhibitors and triazine have been explained by the presence of mutations in target genes, none is responsible for the resistance of this population to glyphosate (Liu et al., 2013a), showing the complexity of NTSR and underlying mechanisms.

The mechanisms of plant responses to chemical stress induced by xenobiotics are not fully understood. Such mechanisms are often studied under conditions of high exposures corresponding to application levels in the field. High levels of xenobiotics strongly impact organism physiology and defense responses,

by inducing molecular injury and damages (Teixeira et al., 2007; Ramel et al., 2009; Das et al., 2010; Nobels et al., 2011; Gomes et al., 2014), mainly related to oxidative stress, membrane disruption, lipid peroxidation, protein damage, or DNA damage. In contrast, few studies deal with the responses to conditions of low exposure to pesticides in a context of runoff or residual contaminations. It has been shown that long-term exposure to sub-lethal pesticide level impacts plant community at the plant development level without causing mortality (Pfleeger et al., 2012). Londo et al. (2014) demonstrated that biomass, flowering phenology, and reproductive functions are affected by sub-lethal glyphosate exposure in *Brassica* spp. Moreover, Ivanov et al. (2013) observed that although sub-lethal concentrations of atrazine did not cause immediate negative and visible effect, long-term exposition impacted the redox homeostasis through an oxidative stress. Long and low herbicide exposure results also in rapid herbicide resistance evolution for exposed populations as demonstrated by Yu et al. (2013) for *Lolium rigidum* in presence of diclofop-methyl. At the molecular level, Das et al. (2010) demonstrated by genome-wide expression profiling that five commercial herbicide formulations at concentration producing a 50% reduction in shoot dry weight (EC50, sub-lethal levels) specifically affected the expression of genes related to ribosome biogenesis and translation, secondary metabolism, cell wall modification and growth. A very recent study demonstrated that subtoxic levels of herbicides acted as chemical hybridization agents, leading to male sterility for the production of hybrid seeds. Their effects were related to reprogramming of gene expression and metabolism in response to low-level herbicide treatments (Li et al., 2015). This study thus showed that complex mechanisms of low-intensity herbicide stress responses may exist. ¹H NMR fingerprinting was also undertaken to analyse substantial metabolic changes in *Lemna minor*'s metabolome after short exposure to different pesticide treatments leading to no phytotoxicity symptom in attempt to develop novel ecotoxicological biomarkers. The discrimination between treatments was mostly based on metabolic variations of substances containing methylene, methine, hydroxy, amine, thiol, olefin, and aldimine groups without mechanistic conclusions (Aliferis et al., 2009). Complex mixtures of xenobiotics at realistic environmental levels produced different effects than those caused by compounds alone or by simple addition of effects (Lydy et al., 2004). Depending on the chemical properties and modes of toxic action of each compound, mixtures can result in higher (synergism) or weaker (antagonism) toxicity, as highlighted for various organisms, as humans, invertebrates or plants (Hertzberg and MacDonell, 2002; Hernández et al., 2012; Cedergreen, 2014). Frankart et al. (2002) have thus shown synergistic effect of a mixture of herbicide (flumioxazin) and copper on photosynthesis inhibition in *L. minor* whereas a mixture of fungicides (fludioxonil or procymidone) and copper produced an antagonism effect. Mixture effects are difficult to analyse and to predict (Dévier et al., 2011; Serra et al., 2013, 2015), and interactions between compounds can alter bioavailability or uptake rate and transport, metabolic activities, target site binding and/or compound excretion (Cedergreen, 2014). Their
study remains however of interest, in particular in the case of no observed effect individual concentrations (Walter et al., 2002).

Hormetic effects and safener effects indicate that xenobiotics can also affect plants under conditions of no adverse effect (NOAE situation: No Observable Adverse Effect) through mechanisms that have seldom been investigated. Hormetic effects that induce beneficial impacts by exposure to low doses of a potentially toxic stressor are achieved through the activation of signal and regulation pathways independently of cellular damage (Velini et al., 2008; Costantini et al., 2010; Belz and Duke, 2014). In that context Nadar et al. (1975) described in Sorghum the growth-promoting effect of atrazine at sub-lethal concentrations in relation with cytokinin-like activity. Stamm et al. (2014) demonstrated in soybean that, even though a thiamethoxam seed treatment did not significantly impacted shoot height and plant biomass, the expression of genes related to plant defense and stress response was altered. Thus, the use of Cruiser® 5FS induces unexpected effects, regarded as cryptic, on a non-target organism. Such cryptic effects were observed in *Arabidopsis thaliana* by Serra et al. (2013) who analyzed the effects of low doses of pesticides, of pesticide degradation products and of their mixtures. In this study, AMPA and hydroxyatrazine, the main degradation products of glyphosate and atrazine, respectively, led to NOAE situations, and nevertheless had significant effects on the expression of genes already known to be affected by high pesticide exposure and on metabolic profiles (Serra et al., 2013). Some chemical treatments induced extensive metabolic changes, such as accumulation of stress-related metabolites (ascorbate) and decrease of carbohydrate levels. Moreover, these chemical stresses effects occurred in parallel with modifications of hormone-related and transcription regulation-related gene expression, thus suggesting underlying regulatory actions of xenobiotic compounds (Serra et al., 2013). Strong interaction effects between chemicals at the molecular, metabolic and physiological levels confirmed that pesticide-related products may act on regulation pathways (Serra et al., 2013).

An another integrative study, focusing on the physiological and metabolic responses of *L. perenne* to diverse subtoxic conditions of chemical stresses, showed primary effects of chemical stressors on seedling metabolism, physiology and growth (Serra et al., 2015). A short exposure to low doses of glyphosate, tebuconazole and their mixture, which consisted of transfer exposure on chemical stressors containing medium during 4 days, did not have any negative effect on root length (NOAE situation), root growth being the most sensitive physiological parameter, but caused cryptic effects on metabolic, regulatory, and signaling processes (Serra et al., 2015). These effects did not however translate into long-term loss of fitness, thus indicating a situation of tolerance to low-level chemical stress. Short exposure was associated with unexpected metabolic changes, as for example significant decrease of Suc and Glc, corresponding to major reorientation of central carbon metabolism. The global analysis in *L. perenne*, combining leaves and root metabolites after short exposure and direct and long exposure, allowed to demonstrate that responses to low chemical stress were associated through a complex

network of metabolic correlations converging on Asn, Leu, Ser, and glucose-6-phosphate (Glc-6-P), which could potentially be modulated by differential dynamics and interconversion of soluble sugars (Suc, trehalose, and Glc; Serra et al., 2015). Such complex metabolic changes reflected chemical stress adjustment rather than deregulation of homeostasis, then leading to root growth maintenance, even under long-term exposure, and thus suggesting the implication of primary mechanisms and molecular regulations (Serra et al., 2015). More importantly, these analyses suggested that complex signaling networks may directly participate in chemical stress responses to rapidly adjust plant metabolism and to counteract mild damaging stress. Moreover, the discovery of major cryptic effects on metabolic, regulatory, and signaling mechanisms under such NOAE conditions raises the issue of the stress concept in plants, as outlined by other authors (Kranner et al., 2010). Within this context, the varying levels of plant stress responses to varying stress intensity may be adaptive. In the absence of harmful effects, a “mild” stress can induce an “alarm response” characterized by post-translational and stress signaling resulting in transcriptomic modifications. The characterization of such low-intensity adaptive mechanisms is therefore of utmost importance, with a number of potential agronomical and ecological applications (stress shield concept).

*L. perenne* is a perennial species of high ecological and agronomic interest. It is one of the predominant forage grasses of high quality in temperate areas (Casler and Duncan, 2003; Barbehenn et al., 2004). Carbon sequestration, soil formation and nutrient cycling are improved by *L. perenne* cover crops (Pouyat et al., 2009). Its value also resides in its relative tolerance to various pollutants of different chemical nature (Dear et al., 2006), which explains why it has been used for phytoremediation strategies (Bidar et al., 2009; D’Orazio et al., 2013).

The present work analyses the molecular responses of *L. perenne*, using a transcriptomic approach, in order to characterize the mechanisms underlying responses to mild chemical stresses induced by no-adverse-effect doses of pesticides. In order to decipher primary mechanisms and pathways involved in adjustments, stress treatments consisted in transfer experiments of non-stressed plants to xenobiotic-containing medium during short exposure. Pesticide treatments consisting of glyphosate, tebuconazole and their combination were applied at NOAE levels to the roots of *L. perenne* seedlings. Both pesticides have been found at residual levels in soils of field margins (Serra et al., 2013) and are frequently detected in runoff water (Potter et al., 2014; Sasal et al., 2015). They are representative examples of agricultural pollution. The impact of root exposure on the whole plant is therefore a primordial process during plant/xenobiotic interactions under conditions of edaphic xenobiotic pollution. Glyphosate is a broad spectrum herbicide (Duke and Powles, 2008) which disrupts the synthesis of aromatic amino acids by inhibiting EPSPS, a key enzyme in the shikimate pathway (Steinrücken and Amrhein, 1980). Tebuconazole, besides acting and being used as a fungicide, can inhibit sterol 14α-demethylase enzymes (Lamb et al., 2001) and limit the rate of gibberellin biosynthesis (Child et al., 1993) in plants. Transcriptomic analysis of *L. perenne* responses to
chemical stress was carried out by RNA-Seq approach involving the pyrosequencing of leaf cDNA libraries (Huang et al., 2012; Ward et al., 2012). A *de novo* transcriptome of *L. perenne* was obtained by the assembly of cDNA reads. A *de novo* assembly was performed using model species data as there is no reference genome available for *L. perenne*. Differential expression of genes and functions were analyzed in parallel with metabolic data in order to obtain a functional insight of molecular regulations induced by low-intensity chemical stress. Our present study shows that, taken together, the effects of pesticide treatments on metabolites and gene expression strongly suggest that low levels of pesticides, whether as single molecule treatment or as mixture, interact, even in the absence of major physiological impact, with plant cell functions as carbohydrate regulations and signaling. The characterization of plant/xenobiotic direct interactions with signaling and hormone cross-talk effects, which is a major field of research in animal toxicology (Frye et al., 2011), should provide novel insights into the environmental impact of low-level runoff or persistent pesticides on plant communities.

# MATERIALS AND METHODS

## Plant Material and Growth Conditions

Seeds of *L. perenne* (Brio cultivar) were briefly washed in ethanol and surface-sterilized in bayrochlore ($20 g L^{-1}$ in water) containing 0.05% tween (v/v) for 20 min and rinsed five times in sterilized water. Moistened seeds were placed in Petri dishes in the dark at 4°C for 7 days in order to break dormancy and homogenize germination. Seeds were sown on pieces of gauze and placed at the top of sterile culture tubes containing liquid growth medium. Gauze pieces were continuously moistened by soaking gauze edges into culture medium, in order to maintain humidity for germination. Germination and growth were carried out under axenic conditions in a control growth chamber at 22°C/20°C under a 16 h light (6000 lux)/8 h dark regime. Growth medium consisted of Hoagland basal salt mix (N°2, Caisson Laboratories, North Logan, UT, USA) adjusted to pH 6. Transfer experiments consisted in xenobiotic exposure of young plants at the same stage of photosynthetic development. After 7 days of growth under control conditions, seedlings were transferred to fresh growth medium containing chemical stressors, thus resulting in xenobiotic exposure at root level. Shoots were harvested 4 days later, corresponding to 11 days of total growth. Different chemical treatments were applied: the broad-spectrum herbicide glyphosate (G, 1 µM), the fungicide tebuconazole (T, 4 µM) and a combination of glyphosate and tebuconazole (GT, 1 µM and 4 µM, respectively). Leaves of seedlings were collected just before the start of the daylight period, ground in liquid nitrogen and stored at -80°C until use.

## Transcriptome Sequencing

Five independent biological replicates, consisting in aerial parts (50 mg fresh weight) of 10 plantlets each, were harvested after transfer experiment and used for transcriptome sequencing. RNA was extracted using TRI Reagent® (Sigma) with an additional DNase treatment. Total RNA samples from each replicate were pooled per treatment condition [Control (C), glyphosate (G),

tebuconazole (T) and glyphosate and tebuconazole mixture (GT)] and resulting samples were polyA-enriched with the Oligotex mRNA kit (Qiagen Cat. N°. 70042). Then, 250 ng of mRNA-enriched RNA were fragmented ($ZnCl_2$), and reverse-transcribed to cDNA using cDNA Synthesis System kit (Roche Cat. N°. 11117831001). For each treatment, one Roche 454 library was prepared and sequenced twice on a quarter of plate on a Roche 454 GS-FLX, using titanium chemistry (titanium chemistry, Sequencing Kit XL+) at the Biogenouest core facility (Rennes, France). The sequencing data are available in the NCBI Sequence Read Archive (SRA) database under the accession reference PRJNA287779 (http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi).

## Assembly and Read Processing

Read filtering and final assembly were performed using the Galaxy instance (Goecks et al., 2010). The 454 reads were processed to remove sequences below 250-bp and above 1000-bp. Reads were improved according to quality and by removing adapter sequences, undetermined bases ("N") and poly-A/T tails. For generation of the reference transcriptome, selected clean reads from all libraries were pooled and assembled to improve mRNA lengths and produce a global library. Reads were assembled in contigs using the *de novo* assembly program Trinity (Haas et al., 2013). Trinity is an Illumina/Solexa-specialized transcriptomic assembler, which is suitable for non-strand-specific and single-end-read data. According to Ren et al. (2012), it gives the best performance among the multiple de Bruijn graph assemblers, and represents an alternative solution for reconstructing full-length transcripts from 454 reads. The assembly was conducted using the default parameters. Reads were ascribed to contigs using RSEM (RNA-seq by expectation-maximization) software and filtered according to the reads per kilobase of target transcript length per million reads mapped (RPKM) values equal to 1 or greater. After a first Trinity assembly, the unmapped (172,634 total unmapped sequences), and unused reads, were submitted again to the Trinity program, leading to a new assembly. Reads that again did not fit into contigs (42,351 unmapped reads) were defined as singletons. These unique sequences were added to the 14,811 unique contigs to constitute the reference transcriptome. The resulting singletons and contigs represented the candidate unigene set.

## Annotation and Functional Analyses

After assembling, tBLASTx alignments (e < 10e-5) against The Arabidopsis Information Resource (TAIR, http://www.arabidopsis.org.gate1.inist.fr) databases were undertaken and unigenes with the highest sequence similarity were functionally annotated using Blast2GO (Conesa et al., 2005) with cut-off e-values of 10e-5 (Blastx) and 10e-6 (mapping). Blast-based annotations were complemented with domain-based annotations using the Inter-ProScan tool (v5). Functional classification of unigenes was based on Blast2GO analysis. Pathway assignments were carried out according to KEGG database (Kyoto Encyclopedia of Genes and Genomes (KEGG) resource; http://www.KEGG.jp/ or http://www.genome.jp/kegg/).
The statistical assessment of Gene Ontology (GO) term enrichments after xenobiotic treatments in comparison to control condition was performed [Fisher's Exact Test with multiple testing correction; false discovery rate (FDR) < 0.05; *p* < 0.005] as implemented in Blast2GO (Conesa et al., 2005). Comparison of conserved amino acids was performed using Clustal Omega for alignment (Sievers et al., 2011).

## Differential Expression Analysis

In order to determine differentially expressed (DE) genes between the 4 conditions, the DESeq method was used. This *in silico* normalization method is included in the DESeq Bioconductor package (Anders et al., 2013). As each condition was represented by a single sample of pooled individuals, the differential expression for each gene was estimated by using the variance of expression for this gene across the 4 conditions. Unigenes were considered as DE at *p* < 0.05.

## qRT-PCR validation

Quantitative RT-PCR was used to confirm *in silico* differential expression and analyse the expression of genes potentially involved in chemical stress response in *L. perenne*. qRT-PCR experiments were carried out using, for each condition, five new independent biological replicates of pooled aerial parts (50 mg fresh weight) from 10 plantlets harvested after transfer experiments. RNA from aerial parts was extracted using TRI Reagent® (Sigma) with an additional DNase treatment. RNA was used for cDNA synthesis (Iscript™ cDNA Synthesis kit, Bio-Rad, Hercules, CA, USA). Resulting cDNAs were used to determine expression profiles according to the different treatments. Quantitative PCR was performed using iQ™ SYBR Green Supermix (Bio-Rad, Hercules, CA, USA). Conditions were as follows: 95°C 3 min, and 40 (95°C 15 s, 60°C 30 s, 72°C 30 s) cycles. All samples were run in duplicate for each primer set. Specific primers for each gene selected for analysis were designed according to 454 sequences using Primer3 software (Rozen and Skaletsky, 2000; Supplemental Table 1). The results of the analysis were treated with Gene Expression version 1.1 software. Relative expressions were assessed in relation to the stable expression level of the GAPDH housekeeping gene (Furet et al., 2012).

## Metabolic Profiling

Five new independent biological replicates of aerial part samples, each consisting in 10 pooled plantlets, were harvested after a transfer experiment, freeze-dried and used for metabolomic profiling. Samples were extracted and analyzed using gas chromatography mass spectrometry (GC/MS) as described by Serra et al. (2015). Metabolite levels were quantified using XCalibur v2.0.7 software (Thermo Fisher Scientific Inc., Waltham, MA, USA) and expressed as nmol.mg⁻¹ of dry weight (DW).

## Statistical Analysis

Metabolic parameters were measured on at least five independent replicates of at least 10 individual plantlets. Gene expression quantification was carried out on at least five other independent biological replicates of 10 pooled plantlets. Statistical analyses

were carried out with version 3.0.1 of R software, analysis between means was carried out using the non-parametric Mann-Whitney-Wilcoxon test.

## RESULTS

### De novo sequence Assembly of the *L. perenne* Shoot Reference Transcriptome

A de novo transcriptome analysis of *L. perenne* was carried out on leaves from plants subjected to a transfer experiment involving NOAE levels of chemical stressors and short periods of exposure. *L. perenne* seedlings were submitted to 3 conditions of short-term low-level xenobiotic exposure [glyphosate (G), tebuconazole (T), and a combination of the two (GT)] in comparison to control condition (C). Leaf cDNA libraries, corresponding to these 4 conditions, were prepared and sequenced using 454 mass sequencing. Before preprocessing, sequencing of cDNA libraries resulted in approximately 0.7 Gbp of sequence data with a GC content of 54%. The genome size of *L. perenne* is estimated to be 2.7 Gb (Fiil et al., 2011). Each quarter plate run was preprocessed, and the runs of each treatment were pooled, thus leading to 300,795, 249,254, 338,011, and 307,033 high-quality reads for, respectively, control, glyphosate, tebuconazole and glyphosate + tebuconazole treatments (Table 1). The reference transcriptome, which was derived from the complete pool of reads from the 4 conditions (C, G, T, GT), contained 1,195,093 reads (Table 1).

Read assembly resulted in an initial reference transcriptome of 14,492 contigs with length ranging from 250 to 6437 bp (Table 2A) and an average contig size of 890 bp. The unmapped reads (172,634 reads) of 250–1000 bp length range were used for a second round of assembly, which resulted in 3564 new contigs. Among the 172,634 reads used, 42,351 remained unmapped and were defined as resulting singletons. Among the 3564 new contigs, 3245 were very similar, if not identical, to contigs built in the first assembly. The 42,351 singletons and the 319 new contigs were added to the 14,492 contigs leading to 57,162 sequences or unigenes (Table 2B) which were used further for whole transcriptome annotation.

The quality of the assembly was checked by comparing the resulting transcriptome to SRA data of *L. perenne* (43,049 sequences, E-MTAB-1556 in ArrayExpress, Vigeland et al., 2013). A tblastx analysis using an *e* < 10e-6 revealed that 45,962 sequences (13,782 contigs and 32,180 singletons) from the present transcriptome matched similar sequences from the deposited SRA data. The reference transcriptome was also compared to the complete mitochondrial (GenBank: JX999996.1) and chloroplastic (NCBI Reference Sequence: NC_009950.1) genomes of *L. perenne* using tblastx program. The mitochondrial genome contains genes corresponding to 14 tRNA, 3 rRNA and 34 proteins (Islam et al., 2013), among which 27 genes matched with high identities (*e* < 10e-6) with unigenes from the present reference transcriptome. Furthermore, 80 genes from the complete chloroplastic genome, which contains genes encoding 76 unique proteins, 30 tRNAs and four rRNAs (Diekmann et al., 2009), were identical or very close (*e* < 10e-6) to the transcriptome unigenes.
<table><caption>TABLE 1 | Summary of 454 sequencing data.</caption><thead><tr><th>Number of high-quality (HQ) reads</th><th></th><th>Control</th><th>Glyphosate</th><th>Tebuconazole</th><th>Glyphosate + Tebuconazole</th></tr></thead><tbody><tr><td rowspan="2">Sequencing reads before preprocessing (Average bp length of HQ reads)</td><td>Run 1</td><td>224,814 (481 ± 184)</td><td>207,483 (450 ± 180)</td><td>234,687 (488 ± 176)</td><td>227,001 (483 ± 183)</td></tr><tr><td>Run 2</td><td>171,316 (401 ± 232)</td><td>143,737 (351 ± 218)</td><td>193,881 (431 ± 224)</td><td>173,328 (411 ± 233)</td></tr><tr><td rowspan="2">Reads after 250–1000 bp trimming</td><td>Run 1</td><td>190,551</td><td>169,605</td><td>202,685</td><td>193,374</td></tr><tr><td>Run 2</td><td>112,834</td><td>84,325</td><td>139,427</td><td>116,918</td></tr><tr><td rowspan="2">Reads after trimming and preprocessing (poly-A/T, N, Quality)</td><td>Run 1</td><td>189,143</td><td>167,578</td><td>200,697</td><td>191,833</td></tr><tr><td>Run 2</td><td>111,652</td><td>81,676</td><td>137,314</td><td>115,200</td></tr><tr><td>Total HQ reads in pools of run 1 and run 2</td><td>Pool</td><td>300,795</td><td>249,254</td><td>338,011</td><td>307,033</td></tr><tr><td>Reference transcriptome</td><td>Pool</td><td></td><td>1,195,093 reads</td><td></td><td></td></tr></tbody></table>

<table><caption>TABLE 2 | Assembly results.</caption><thead><tr><th></th><th>Reference</th><th>Control</th><th>Glyphosate</th><th>Tebuconazole</th><th>Glyphosate + Tebuconazole</th></tr></thead><tbody><tr><td colspan="6"><strong>(A) FIRST ASSEMBLY</strong></td></tr><tr><td>Number of contigs</td><td>14,492</td><td>10,155</td><td>5773</td><td>9109</td><td>9408</td></tr><tr><td>Total length of contigs (bp)</td><td>12,905,107</td><td>10,070,689</td><td>6,098,760</td><td>9,166,250</td><td>9,611,071</td></tr><tr><td>Number of contig ≥ 1000 bp</td><td>4342</td><td>3937</td><td>2549</td><td>3661</td><td>3921</td></tr><tr><td>Mean contig length (bp)</td><td>890 ± 466</td><td>991 ± 498</td><td>1056 ± 537</td><td>1006 ± 508</td><td>1031 ± 503</td></tr><tr><td colspan="6"><strong>(B) FINAL ASSEMBLY</strong></td></tr><tr><td>First assembly</td><td>14,492</td><td>8089</td><td>3646</td><td>7005</td><td>7256</td></tr><tr><td>Second assembly</td><td>319</td><td>130</td><td>50</td><td>94</td><td>103</td></tr><tr><td>Singletons</td><td>42,351</td><td>8089</td><td>6218</td><td>12,553</td><td>11,025</td></tr><tr><td>Total</td><td>57,162</td><td>22,087</td><td>9914</td><td>19,652</td><td>18,384</td></tr></tbody></table>

In order to validate further the assembled unigenes, sequence-based alignments were performed against the TAIR (Arabidopsis thaliana) database by using the tblastx algorithm (Altschul et al., 1997). According to tblastx data, 55.44% (31,696 unigenes) of the matched sequences showed strong homology with TAIR data ($e < 10e-50$), and 84.9% (48,550 unigenes) of the top hits exhibited lower but significant homology ($e$-value over $10e-6$).

## Functional Annotation of the *L. perenne* Shoot Reference Transcriptome

Functional annotation of the reference transcriptome was undertaken using the Blast2GO tool and TAIR database. 9311 (16% of all unigenes) unigenes showed no blast hit in the TAIR (Arabidopsis thaliana) database. Among the unigenes yielding blast results, 204 (0.32%) were not mapped and 1296 (2.2%) mapped unigenes were not annotated, thus leading to 46,330 (81%) annotated unigenes. In many cases, multiple GO terms were assigned to the same unigene.

GO terms were classified into functional groups according to biological process, molecular function, and cellular component classes. Distribution of the various GO terms identified for the *L. perenne* reference transcriptome is presented in **Figure 1** for the major GO classes. The proportions of biological activities were similar to those described in other *L. perenne* transcriptomes (Farrell et al., 2014; Duhoux et al., 2015). GO terms linked to metabolism and biosynthetic processes represented more than 50% of the biological processes (**Figure 1A**). GO terms related to stress responses, signaling and regulations (response to stress, response to external/abiotic/chemical stimulus, detection of stimulus, single organism signaling, regulation of biological processes) were also well represented (16.6%). Concerning molecular functions (**Figure 1B**), GO terms linked to various binding activities, as well as activities related to translation (structural constituent of ribosomes), were highly represented. Functional categories also covered enzyme activities such as hydrolase, transferase, oxidoreductase, or lyase activities. The subcellular localization of these processes and functions (**Figure 1C**) showed the importance of gene products associated with intracellular organelles, chloroplasts (thylakoids, light-harvesting complexes), cytoplasm and ribosomal complexes. Functional distribution therefore demonstrated that this *L. perenne* reference transcriptome exhibited a good coverage of essential types of plant biological activities.
The relative importance of various metabolic pathways was assessed from the assignment of unigene annotations in the KEGG database (Figure 2). Most of the predominant pathways were found to belong to carbon metabolism (Glyoxylate and dicarboxylate metabolism, Glycolysis/Gluconeogenesis, Starch and sucrose metabolism, Pentose phosphate pathway, Fructose and mannose metabolism, Amino sugar and nucleotide sugar metabolism, Pyruvate metabolism, TCA cycle, Galactose metabolism, Ascorbate, and aldarate metabolism) and to energy metabolism (Carbon fixation in photosynthetic organisms, Methane metabolism, Oxidative phosphorylation, Nitrogen metabolism, Photosynthesis). Five pathways related to amino acid metabolism were also represented.

## Effects of Short-Term Root-Level Exposure to NOAE Levels of Chemical Stressors on the Relative Regulation of Functional Categories in *L. perenne* leaves

### Assessment of Xenobiotic-Related Molecular Effects

The relative proportions of functional categories, as defined by GO terms, were analyzed relative to the different conditions of chemical stress. GO enrichment analysis for each class (biological processes, cellular components, and molecular functions) revealed the GO terms for which the quantity of annotated unigenes (contigs and singletons) presented significant differences (FDR < 0.05 and *p* < 0.005) after xenobiotic treatments in comparison to control condition. Results were expressed as the percentage of unigenes in each GO term category in a given condition relative to the total number of annotated unigenes in the same class. The distribution of significantly-enriched GO term categories is given in Figure 3. The percentage of annotated reads in each GO term category in a given condition was also calculated relative to the total number of annotated reads in the same class. In this analysis of annotated read levels, only GO term categories for which at least 50 reads had been counted in at least one condition were selected (Supplemental Table 2). Changes in proportions of GO term categories in terms of unigene enrichment were likely to reflect the effects of xenobiotic treatments on global regulatory processes of genome expression affecting each category. Modifications of proportions of GO term categories in terms of read proportion gave further information on induction or repression.

In a number of cases, the relative proportions of GO term categories were decreased by xenobiotic treatments in terms of percentages relative to annotated unigenes (Figure 3), and also in terms of percentages relatively to annotated reads (Supplemental Table 2). Nevertheless, among all the significant differences induced by treatments, G and T treatments mainly negatively affected the proportion of GO terms relative to control. By contrast, GT treatment induced as much enrichment as depletion in the proportions of GO terms relative to control (Figure 3). Such differences strongly suggested that the mixed GT treatment induced specific effects in comparison with those of single xenobiotic treatments as “response to monosaccharide stimulus”, “sugar mediated signaling pathway”, “starch catabolic process”, “negative regulation of translation”, “pentose phosphate shunt”,

A

<figure><img src="image_3.png" /></figure>

Biological
processes

B

<figure><img src="image_4.png" /></figure>

Molecular
functions

C

<figure><img src="image_5.png" /></figure>

Cellular
components

**FIGURE 1 | Distribution of Gene Ontology (GO) terms in the reference shoot transcriptome of *L. perenne*.** Biological processes (A), Molecular functions (B) and Cellular components (C) classes are shown. The relative size of pie chart sectors reflects the percentage of annotated unigenes of each GO term category.
<figure><img src="image_6.png" /><figcaption>**FIGURE 2** | Distribution of metabolic pathways in the reference shoot transcriptome of *L. perenne*. Assignment of annotated unigenes to metabolic pathways was carried out with the KEGG database (Kyoto Encyclopedia of Genes and Genomes (KEGG) resource; http://www.kegg.jp/ or http://www.genome.jp/kegg/). Results are expressed as percentages of annotated unigenes of each metabolic pathway.</figcaption></figure>

"pyrophosphatase activity", "NADP metabolic process", "ATP binding", "Golgi apparatus part". Moreover G treatment often led to more pronounced effects than those of T and GT treatments (Figure 3).

## Effects on Gene Expression and Protein Translation Processes

The "structural constituent of ribosome" category (Figure 3B) showed increased proportions in the presence of G. The potential involvement of ribosome biogenesis in xenobiotic responses was also reflected by the significant increase of "ribosomal subunit" category in the presence of xenobiotics and particularly in presence of G (Figure 3C). Other aspects of the biogenesis of functional ribosomes seemed to be sensitive to chemical stresses, since categories of linked GO terms such as "ribonucleoprotein complex biogenesis", "rRNA processing" and "ribosome biogenesis" were negatively affected by chemical stress treatments. G and GT treatments also induced a slight decrease of unigene proportion for "DNA-dependent transcription" (Figure 3A). The stability of unigene proportion related to "negative regulation of translation" in the presence of T and G contrasted with a strong decrease of unigene proportion in presence of GT, thus suggesting combination-specific regulations (Figure 3A). For upper-level biological process classes, such as "translation", overrepresented in presence of G and T, and "gene expression", overrepresented in presence of G and slightly

underrepresented in presence of GT, the increase of unigene proportion in the presence of G (Figure 3A) was associated with an increase in read number (Supplemental Table 2), thus indicating the induction of specific genes. This increase of "gene expression" category occurred in parallel with a decrease in the "gene silencing" category (Figure 3A). Finally, negative effects of G on "histone modification"-related unigenes were observed in upper levels of classification such as "chromatin organization" and "macromolecule modification" (Figure 3A).

## Effects on Photosynthesis and ATP Dynamics Processes

Unigenes linked to the "ribulose-1,5-bisphosphate carboxylase/oxygenase activator activity" or "Rubisco activator" showed decreased proportions in response to xenobiotics, particularly in response to G (Figure 3B). The associated "ribulose-bisphosphate carboxylase activity" category was significantly more represented in response to T and underrepresented in the presence of G and GT. G and GT treatments reduced unigene proportions corresponding to biological processes associated with "carbon fixation", "photosynthesis", "photosystem II assembly", and "photosynthetic electron transport chain" (Figure 3A). Glyphosate application also reduced the number of unigenes involved in "chlorophyll biosynthetic process" and in related pathways such as "porphyrin-containing compound metabolic
**A** Biological processes

<figure><img src="image_7.png" /></figure>

**B** Molecular functions

<figure><img src="image_8.png" /></figure>

**C** Cellular components

**FIGURE 3** | Effects of xenobiotics on the distribution of representatively-enriched Gene Ontology (GO) terms in the shoot transcriptomes of *L. perenne*. Biological processes (A), Molecular functions (B), and Cellular components (C) classes are shown. Results are expressed as percentages of enriched and annotated unigenes for each GO term category and for each condition (control, glyphosate, tebuconazole, glyphosate plus tebuconazole) relatively to the total number of annotated unigenes in that class. Stars indicated FDR < 0.05 and *p* < 0.005 after Fisher's Exact Test with robust FDR (false discovery rate) correction.
process” and “tetrapyrrole biosynthetic process” (Figure 3A). While no unigene enrichment of “magnesium chelatase activity” category was observed, there was an increase in read numbers in the presence of the 3 xenobiotic treatments, for genes whose functions were related to these chlorophyll biosynthetic pathway enzymes (Supplemental Table 2). Xenobiotic-related perturbations in chloroplastic and photosynthetic pathways were also identified for cellular component GO terms (Figure 3C). Numerous unigenes were ascribed to GO terms related to “chloroplast”, “thylakoid”, “plastoglobule”, “chloroplast thylakoid”, and “thylakoid membrane”. These categories generally showed higher proportion of unigenes in the presence of T and GT, while these xenobiotics tended to decrease slightly the proportion of unigenes related to “cytochrome b6f complex” categories (Figure 3C).

The “binding” GO term category was well represented (Figure 3B), in particular with regard to “nucleoside phosphate binding”, “ATP binding”, and “ADP binding”. Proportions of annotated unigenes for these molecular functions decreased in the presence of G, and slightly increased in the presence of GT. However, concerning the “ATP binding” annotation, which was highly represented among GO terms, an increase of the number of corresponding reads was observed in the presence of G, suggesting strong expression of specific and constitutive unigenes (Supplemental Table 2). Likewise, other functions related to nucleoside and nucleotide dynamics, which showed frequent occurrences and high percentages, were responsive to xenobiotic treatments. Thus, the GO term category “nucleoside-triphosphatase activity” showed decreased proportion in the presence of G, whereas GT treatment increased it (Figure 3B). This metabolic activity was characterized by an increase of related read numbers in the presence of G, T, and GT, suggesting strong expression of specific unigenes (Supplemental Table 2). This pattern of modifications was reflected in the changes affecting the corresponding upper-level GO term, “hydrolase activity”. Nucleoside-triphosphatase activity is an essential provider of energy for active transport and is related to the “P-P-bond-hydrolysis-driven transmembrane transporter activity” category, which showed higher proportions of unigenes in the presence of GT (Figure 3B).

## Effects on Environmental Response Processes

G and GT treatments induced slight decreases of unigene proportions for many biological processes related to abiotic and biotic stress signaling processes, such as “jasmonic acid mediated signaling pathway”, “regulation of plant-type hypersensitive response”, “regulation of programmed cell death”, and “cellular response to stress” (Figure 3A). Moreover, the proportions of unigenes and reads related to “response to salicylic acid”, “MAPK cascade”, “intracellular protein kinase cascade”, and “regulation of hydrogen peroxide metabolic process” were negatively affected by all the chemical stressors (Figure 3A). Most of the biological processes linked to light responses showed a decrease for the proportions of related unigenes in the presence of G and GT. GO term categories such as “response to light stimulus”, and related sublevels such as “non-photochemical quenching”, “response to blue light”, “response to red or far red light”, and “response to high light intensity”, followed the same tendency (Figure 3A).

## Effects on Metabolic Functions

The proportions of unigene categories linked to carbohydrates and their derivatives were significantly influenced by chemical treatments. Unigenes annotated in “starch catabolic process”, “sugar mediated signaling pathway”, and “response to monosaccharide stimulus” in biological processes were induced by GT treatments (Figure 3A). GT condition reduced the number of unigenes related to “pentose-phosphate shunt” (Figure 3A). The effects of chemical treatments on the “pentose-phosphate shunt” category were reflected in the responses of closely related categories such as “NADPH regeneration”, “NADP metabolic process”, “nicotinamide nucleotide metabolic process”, “pyridine nucleotide metabolic process”, “oxidoreduction coenzyme metabolic process”, “coenzyme metabolic process”, and “cofactor metabolic process” (Figure 3A).

A number of transport activities in the molecular function category were found to be potentially induced by xenobiotics, particularly by T. This was the case for the “organophosphate ester transmembrane transporter” and “carboxylic acid transmembrane transporter activity” categories (Figure 3B). There was also an increase in the unigene proportion of “cytoplasmic membrane-bound vesicle” category in response to GT (Figure 3C).

Finally, GT treatment decreased the unigene proportion of “delta1-pyrroline-5-carboxylate synthetase (P5CS) activity” (Figure 3B) implicated in the proline biosynthesis pathway.

## Effects on Growth and Development Processes

All chemical treatments, and more particularly T and G, increased the number of unigenes annotated in the “adenosylmethionine decarboxylase (SAMDC) activity” molecular function category (Figure 3B), which is related to spermidine and spermine biosynthesis from putrescine. In contrast, for the “polyamine oxidase activity” category, which catalyzes the oxidative degradation of polyamines, the related unigenes were not differentially enriched while the number of reads decreased in response to T, G, and GT (Supplemental Table 2). Such maintenance of polyamine levels, which are essential for growth and development in plants (Galston and Sawhney, 1990), may be, at least in part, related to the variations of categories at a higher level of organization. Indeed, T and GT treatments enhanced the number of unigenes correlated to developmental processes such as “seed maturation”, “primary shoot apical meristem specification” and more generally “regulation of cell differentiation”, while G treatment had opposite effects on proportions of unigenes and reads (Figure 3A, Supplemental Table 2). The presence of G also induced a decrease of unigene and read proportions linked to “leaf development” and “organ morphogenesis”.

## Identification and Characterization of Xenobiotic-Responsive Genes in *L. perenne*

Using a *p* < 0.05, DESeq analysis revealed that 69 unigenes were significantly DE in at least one of the four conditions (Table 3). These gene expression profiling data obtained from
<table><caption>TABLE 3 | Significantly differentially expressed transcripts (DESeq analysis).</caption><thead><tr><th rowspan="2">Unigene ID</th><th rowspan="2">Description (TAIR)</th><th colspan="3">DEseq analysis Log<sub>2</sub> [ratio]</th></tr><tr><th>T/C</th><th>G/C</th><th>GT/C</th></tr></thead><tbody><tr><td>comp7278_c0_seq1</td><td>Magnesium-chelatase subunit chlH</td><td>Repressed</td><td>Repressed</td><td>0.304</td></tr><tr><td>comp7493_c0_seq1</td><td>Protein FLUORESCENT IN BLUE LIGHT (FLU)</td><td>-2.456</td><td>Repressed</td><td>-1.124</td></tr><tr><td>comp7590_c0_seq3</td><td>Plastid transcriptionally active (PTAC16)</td><td>Induced</td><td>Induced</td><td>Induced</td></tr><tr><td>comp7590_c0_seq4</td><td>Plastid transcriptionally active (PTAC16)</td><td>Repressed</td><td>Repressed</td><td>-1.124</td></tr><tr><td>comp7427_c0_seq4</td><td>Light harvesting complex photosystem II</td><td>Repressed</td><td>Repressed</td><td>-1.149</td></tr><tr><td>comp7634_c0_seq2</td><td>Photosystem I, PsaA/PsaB protein</td><td>0.939</td><td>Repressed</td><td>0.069</td></tr><tr><td>comp7693_c0_seq11</td><td>Photosystem II reaction center protein D</td><td>3.129</td><td>4.874</td><td>Repressed</td></tr><tr><td>comp7726_c0_seq27</td><td>Photosystem II reaction center protein B</td><td>-1.363</td><td>Repressed</td><td>-2.879</td></tr><tr><td>comp7726_c0_seq5</td><td>Photosystem II reaction center protein B</td><td>-1.704</td><td>1.167</td><td>Repressed</td></tr><tr><td>comp7726_c0_seq1</td><td>Photosystem II reaction center protein B</td><td>-0.434</td><td>Repressed</td><td>-0.654</td></tr><tr><td>comp6560_c0_seq1</td><td>Photosystem I P subunit</td><td>1.222</td><td>Repressed</td><td>-1.372</td></tr><tr><td>comp6663_c0_seq5</td><td>Isoleucyl-tRNA synthetase-like protein</td><td>Repressed</td><td>Repressed</td><td>1.706</td></tr><tr><td>comp7674_c0_seq3</td><td>Phosphoenolpyruvate carboxylase</td><td>Induced</td><td>NDE</td><td>NDE</td></tr><tr><td>comp6943_c0_seq1</td><td>Uncharacterized protein</td><td>-0.380</td><td>Repressed</td><td>-0.366</td></tr><tr><td>comp6984_c0_seq1</td><td>Rubredoxin family protein</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp6746_c0_seq2</td><td>Actin 7</td><td>Repressed</td><td>2.586</td><td>3.379</td></tr><tr><td>comp7483_c1_seq2</td><td>Phosphoglycerate kinase family protein</td><td>-0.723</td><td>Repressed</td><td>-2.239</td></tr><tr><td>comp7355_c0_seq2</td><td>CIPK9, SnRK3.12, PKS6, CBL-interacting protein kinase</td><td>1.097</td><td>Repressed</td><td>1.794</td></tr><tr><td>comp7489_c0_seq2</td><td>Calcium sensing receptor</td><td>Repressed</td><td>2.415</td><td>0.654</td></tr><tr><td>comp7362_c0_seq6</td><td>Aldehyde dehydrogenase 3F1</td><td>-2.379</td><td>2.737</td><td>0.291</td></tr><tr><td>comp7713_c0_seq4</td><td>Protein kinase superfamily protein</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7705_c0_seq3</td><td>Phototropin serine-threonine kinase</td><td>Repressed</td><td>Repressed</td><td>0.645</td></tr><tr><td>comp7527_c0_seq1</td><td>Protein phosphatase 2C family protein</td><td>2.737</td><td>Repressed</td><td>-0.516</td></tr><tr><td>comp7476_c0_seq2</td><td>Ferredoxin-NADP(+)-oxidoreductase 1 (FNR)</td><td>-1.100</td><td>Repressed</td><td>0.670</td></tr><tr><td>comp7633_c0_seq4</td><td>4-hydroxy-3-methylbut-2-enyl diphosphate synthase</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7402_c0_seq2</td><td>F-type H+-transporting ATPase subunit delta</td><td>Repressed</td><td>4.265</td><td>1.863</td></tr><tr><td>comp7682_c1_seq3</td><td>H(+)-ATPase 2</td><td>Repressed</td><td>-1.322</td><td>0.192</td></tr><tr><td>comp7721_c0_seq19</td><td>ATP synthase subunit alpha</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7721_c0_seq17</td><td>ATP synthase subunit alpha</td><td>Repressed</td><td>1.415</td><td>-1.294</td></tr><tr><td>comp7700_c0_seq2</td><td>ABC transporter F family</td><td>-0.214</td><td>Repressed</td><td>0.016</td></tr><tr><td>comp7571_c0_seq5</td><td>Major facilitator superfamily protein</td><td>0.959</td><td>Repressed</td><td>0.564</td></tr><tr><td>comp7726_c0_seq35</td><td>RNA polymerase subunit alpha</td><td>Induced</td><td>NDE</td><td>NDE</td></tr><tr><td>comp7697_c0_seq14</td><td>Homeodomain-like superfamily protein</td><td>-0.029</td><td>Repressed</td><td>0.232</td></tr><tr><td>comp7697_c0_seq11</td><td>Homeodomain-like superfamily protein</td><td>0.222</td><td>Repressed</td><td>Repressed</td></tr><tr><td>comp7688_c0_seq2</td><td>NAD(P)-binding Rossmann-fold superfamily protein</td><td>Repressed</td><td>0.023</td><td>-0.101</td></tr><tr><td>comp7592_c0_seq2</td><td>Translation initiation factor IF</td><td>-0.778</td><td>Repressed</td><td>1.128</td></tr><tr><td>comp6892_c0_seq2</td><td>Gametogenesis related family protein</td><td>0.824</td><td>2.314</td><td>0.914</td></tr><tr><td>comp7302_c0_seq1</td><td>Uncharacterized protein</td><td>0.740</td><td>1.816</td><td>0.545</td></tr><tr><td>comp7589_c0_seq4</td><td>Chloroplast stem-loop binding protein</td><td>3.392</td><td>4.265</td><td>Repressed</td></tr><tr><td>comp7452_c0_seq3</td><td>Indoleacetic acid-induced protein</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7604_c0_seq2</td><td>Tubulin alpha-5</td><td>-2.100</td><td>Repressed</td><td>0.429</td></tr><tr><td>comp7501_c0_seq5</td><td>DNA photolyase</td><td>-0.314</td><td>Repressed</td><td>-0.982</td></tr><tr><td>comp7711_c0_seq1</td><td>Presequence protease 1</td><td>-1.397</td><td>-2.426</td><td>Repressed</td></tr><tr><td>comp7711_c0_seq5</td><td>Presequence protease 1</td><td>NDE</td><td>Induced</td><td>Induced</td></tr><tr><td>comp7711_c0_seq3</td><td>Presequence protease 1</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7458_c0_seq2</td><td>Peroxidase superfamily protein</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7458_c0_seq1</td><td>Peroxidase superfamily protein</td><td>Repressed</td><td>-0.285</td><td>-1.409</td></tr><tr><td>comp7628_c1_seq8</td><td>Senescence-associated protein</td><td>-0.555</td><td>1.415</td><td>Repressed</td></tr><tr><td>comp7309_c0_seq3</td><td>Glycolate oxidase 1 (GOX1)</td><td>2.410</td><td>Repressed</td><td>-1.294</td></tr><tr><td>comp7324_c0_seq3</td><td>Phenylalanine ammonia-lyase</td><td>Repressed</td><td>0.337</td><td>-0.957</td></tr></tbody></table>

(Continued)
<table><caption>TABLE 3 | Continued</caption><thead><tr><th rowspan="2">Unigene ID</th><th rowspan="2">Description (TAIR)</th><th colspan="3">DEseq analysis Log<sub>2</sub> [ratio]</th></tr><tr><th>T/C</th><th>G/C</th><th>GT/C</th></tr></thead><tbody><tr><td>comp7722_c0_seq7</td><td>Phenylalanine ammonia-lyase</td><td>Repressed</td><td>1.678</td><td>Repressed</td></tr><tr><td>comp7722_c0_seq5</td><td>Phenylalanine ammonia-lyase</td><td>-0.710</td><td>-2.977</td><td>Repressed</td></tr><tr><td>comp7722_c0_seq3</td><td>Phenylalanine ammonia-lyase</td><td>NDE</td><td>NDE</td><td>Induced</td></tr><tr><td>comp6929_c0_seq4</td><td>O-methyltransferase</td><td>Repressed</td><td>0.093</td><td>-1.031</td></tr><tr><td>comp7591_c1_seq3</td><td>ACC oxidase 5 (ACO5)</td><td>Induced</td><td>Induced</td><td>Induced</td></tr><tr><td>comp7591_c1_seq4</td><td>ACC oxidase 5 (ACO5)</td><td>Repressed</td><td>-0.700</td><td>Repressed</td></tr><tr><td>comp6832_c0_seq2</td><td>β-1,3-endoglucanase</td><td>2.485</td><td>Repressed</td><td>0.969</td></tr><tr><td>comp7382_c0_seq3</td><td>Glutamate synthase 1 [NADH]</td><td>0.753</td><td>Repressed</td><td>-0.294</td></tr><tr><td>comp7386_c0_seq10</td><td>Polyamine oxidase</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7536_c0_seq5</td><td>Sulfite reductase</td><td>Induced</td><td>NDE</td><td>Induced</td></tr><tr><td>comp7497_c0_seq1</td><td>ATP sulfurylase 2</td><td>Repressed</td><td>1.830</td><td>-0.071</td></tr><tr><td>comp7497_c0_seq2</td><td>ATP sulfurylase 2</td><td>0.222</td><td>Repressed</td><td>Repressed</td></tr><tr><td>comp6892_c0_seq1</td><td>Uncharacterized protein</td><td>0.850</td><td>2.148</td><td>0.842</td></tr><tr><td>comp7454_c0_seq1</td><td>Uncharacterized protein</td><td>-0.915</td><td>1.126</td><td>Repressed</td></tr><tr><td>comp7706_c0_seq8</td><td>Dentin sialophosphoprotein-related</td><td>Repressed</td><td>0.699</td><td>Repressed</td></tr><tr><td>comp7528_c0_seq1</td><td>Uncharacterized protein</td><td>2.170</td><td>1.608</td><td>-0.346</td></tr><tr><td>comp7706_c0_seq6</td><td>Dentin sialophosphoprotein-related</td><td>Induced</td><td>Induced</td><td>Induced</td></tr><tr><td>comp7188_c0_seq3</td><td>Unknown protein</td><td>Induced</td><td>Induced</td><td>Induced</td></tr><tr><td>comp6657_c0_seq2</td><td>Pentapeptide repeat-containing protein</td><td>Induced</td><td>Induced</td><td>Induced</td></tr></tbody></table>

C, control; G, glyphosate; T, tebuconazole; GT, glyphosate + tebuconazole. The log<sub>2</sub> [ratio] values could not be calculated when no expression was detected for a given transcript. Such cases were indicated by "Induced" when no expression was detected in the control, "Repressed" when no expression was detected in the treatment conditions. Unigenes were considered as differentially expressed at *p* < 0.05. "NDE" indicated cases for which transcripts were non-differentially expressed. Unigene descriptions were established according to best blast hits in TAIR. The intensities of red or green colors respectively indicate increase or decrease of expression levels in the presence of G, T, and GT compared to the control condition.

RNA-seq were complemented with qRT-PCR analysis of the relative expression of candidate genes (Figure 4). The alignments of corresponding protein sequences of candidate unigenes with *Arabidopsis thaliana*, *Oryza sativa*, and *Brachypodium distachyon* orthologs (Supplemental Figure 1) highlighted high degree of sequence conservation between *L. perenne* and these monocot and dicot model species, and strengthened the automatic *in silico* annotation.

Most of the observed differential expressions (Table 3) validated the analysis of GO term category variations under conditions of xenobiotic treatment (Figure 3). Among the DE unigenes, a high proportion of xenobiotic-repressed genes were associated with photosynthesis or chlorophyll biosynthetic processes. This was particularly the case for the transcripts of genes related to the tetrapyrrole pathway, such as those encoding Magnesium-chelatase subunit chlH and Fluorescent In Blue Light (FLU), a tetratricopeptide repeat (TPR)-containing protein, highly conserved between *A. thaliana*, *O. sativa*, and *B. distachyon* (Supplemental Figure 1) that were decreased by xenobiotic treatments (Table 3). This decrease was confirmed by qRT-PCR, when glyphosate or tebuconazole were applied alone (Figure 4). Similar expression profiles were observed for all photosynthesis- and chlorophyll-related unigenes, except for one isoform encoding the plastid transcriptionally active 16 (PTAC16) gene that was induced by all conditions (Table 3). Induction of expression under G and T treatments was detected

by DESeq for the transcript annotated "photosystem II reaction centre protein D". This was also the case for some genes related to photosystems, which exhibited significant induction under either G or T treatments (Table 3).

The DESeq analysis revealed that a high proportion of xenobiotic-responsive genes are involved in signaling, and more particularly, correspond to signal transduction through protein-kinases and protein-phosphatases (Table 3). This was the case for the CIPK9-annotated gene, which encodes a Calcineurin B-Like Protein (CBL)-interacting protein kinase highly conserved between monocot and dicot species (Supplemental Figure 1). This Serine/Threonine protein kinase annotated unigene, also known as sucrose non-fermenting (SNF)-related kinase family of serine/threonine kinase 3.12 (SNRK3.12), as well as a calcium sensing receptor-annotated transcript, were found to be DE (Table 3). DESeq showed a repression of CIPK9 by G and an induction by T and GT (Table 3). qRT-PCR analysis confirmed this repression by G (Figure 4).

Transcripts related to protein phosphorylation and dephosphorylation were highly represented in the different datasets. The effects of xenobiotic stresses on highly conserved Serine/Threonine protein kinases (Supplemental Figure 1), such as comp7296_c0_seq6, annotated as a Serine/Threonine protein kinase, and a unigene comp2769_c0_seq1, annotated as a SNRK2 protein (Kertesz et al., 2002), were analyzed by
<figure><img src="image_9.png" /><figcaption>**FIGURE 4** | Effects of xenobiotics on the transcript levels of candidate xenobiotic-responsive genes in leaves of *L. perenne*. *L. perenne* plants were grown for 7 days in the absence of pesticide and then transferred to fresh growth medium containing the different xenobiotic treatments [glyphosate (G, 1 µM), tebuconazole (T, 4 µM) and a combination of glyphosate and tebuconazole (GT, 1 µM and 4 µM, respectively)]. Transcript levels (mean ± SEM) were quantified by qRT-PCR after 4 days of treatment. Unigene descriptions were established according to best blast hits in TAIR. Statistical analysis between means was carried out using the Mann-Whitney-Wilcoxon test. Statistical significance of differences (*p* < 0.05) between treatments is indicated by different letters above bars.</figcaption></figure>

qRT-PCR (Figure 4). Glyphosate, like for the CIPK9-annotated gene, repressed expression of these genes (Figure 4). As observed for SnRK2-type unigenes, glyphosate tended to repress a PP2C-annotated gene, comp7527_c0_seq1 (Table 3, Figure 4). Different types of expression profile were found, such as the high upregulation of F-type H⁺-transporting ATPase subunit delta-encoding gene by G, the specific repression of the unigene related to H⁺-ATPase 2 by G or T, or the induction by T or GT of one of the isoforms encoding an ATP synthase subunit alpha (Table 3).

The expression of unigenes associated with the regulation of translation or transcription also differed depending on the chemical stressor. Glyphosate repressed the expression of genes associated with Homeodomain-like superfamily proteins and involved in the regulation of transcription (Table 3). Another transcription elongation factor encoded by comp5339_c0_seq1 presented the same kind of expression profile (Figure 4). On the contrary, the expression of 3 unigenes, whose function was linked to translation (comp6892_c0_seq2, comp7302_c0_seq1, and comp7589_c0_seq4), was significantly
upregulated by glyphosate and GT (Table 3). A GOX1 gene encoding a Glycolate oxidase was differentially repressed in the presence of glyphosate and induced in presence of tebuconazole (Table 3). Genes encoding Cytochrome P450 reductase proteins tended to be repressed by glyphosate, whereas the GT mixture had opposite effects (Figure 4).

The qRT-PCR analysis of one of genes annotated as Phenylalanine ammonia-lyases (PAL) demonstrated the repressive effect of xenobiotics on its expression (Figure 4). A gene homologous to a gene encoding a caffeate O-methyltransferase was repressed by T and GT (Table 3). 1-aminocyclopropane-1-carboxylic acid (ACC) oxidases were characterized in the shoot transcriptome of *L. perenne* by two isoforms that were differentially regulated, with one being induced by chemical stresses, while the other was repressed (Table 3). In parallel, the gene annotated as β-1,3-endoglucanase was induced by T and GT, as confirmed by qRT-PCR trends for T treatment (Table 3, Figure 4). The sulfur assimilation pathway, which leads to cysteine biosynthesis and *in fine* to glutathione (GSH) biosynthesis, involves *inter alia*, ATP sulfurylase and sulfite reductase (Anjum et al., 2015). In *L. perenne*, one of the 2 isoforms annotated as ATP sulfurylase showed repression by T and induction by G, while the other was repressed by G and GT (Table 3). The qRT-PCR analysis also showed that G tended to repress this latter isoform (Figure 4). In contrast, T and GT had a positive impact on the transcription of the sulfite reductase-annotated gene (Table 3). Finally, differential expression analysis revealed several uncharacterized transcripts whose expression, for most of them, was induced in the presence of glyphosate and with contrasting effects of the other treatments (Table 3).

## Metabolic Modifications in *L. perenne* leaves Under Conditions of Short-Term Root-Level Exposure to NOAE Levels of Chemical Stressors

The effects of the present xenobiotic stress conditions (low intensity, short-term exposure, root exposure) were also analyzed in terms of metabolic changes in *L. perenne* leaves (Figure 5). Moreover, as previously described in *A. thaliana* (Serra et al., 2013), major metabolic changes were found to occur despite the absence of major stress perturbations. All of the subtoxic conditions analyzed here (G, T, GT) induced or tended to induce a decrease of soluble sugar contents (Fru, Glc, Suc) whereas Tre levels were no affected (Figure 5). In contrast, G and T acted differently on the levels of Fru-6-P compared to the GT mixture, with an increase of this phosphorylated sugar with G and T and no effect with GT. Glc-6-P remained relatively stable in response to the three conditions.

Some metabolic changes could be related to the patterns of differential expression described above (Figures 3, 4, Table 3). Glyphosate or tebuconazole exposure induced a decrease of glycine levels (Figure 5) and also affected the expression of glycine-related genes such as glycine decarboxylase (Figure 4). The decrease in putrescine levels in *L. perenne* leaves in response

to xenobiotic treatments (Figure 5) was in line with the induction by T and GT of the unigene encoding a polyamine oxidase (Table 3) and with significant variations of related GO term categories, as “adenosylmethionine decarboxylase (SAMDC) activity” and “spermine biosynthetic process” (Figure 3). The decrease of glutamate levels in *L. perenne* leaves in response to G (Figure 5) was reminiscent of the decrease of glutamate synthase expression (Table 3).

All treatments tended to increase the levels of arabinose, Phe, Lys, Leu, Ile, Val, Tyr, citrate and lactate and to decrease the levels of putrescine and Pro. In particular, glyphosate treatment did not decrease levels of aromatic amino acids (Trp, Tyr, Phe). For Ala, Asp, Fru-6-P, glutamate and ornithine, effects of G or T alone were mostly comparable, and GT mixture effects were not additive (Figure 5). In the case of stress metabolites, pipecolate, putrescine, inositol and Pro decreased or were not affected, while lactate increased by low chemical stress (Figure 5). This decrease of Pro levels of plants subjected to GT (Figure 5) could be related to the reduction of proportion of genes related to “P5CS activity” category (Figure 3B).

# DISCUSSION

## Xenobiotic-Regulated RNA-Seq and qRT-PCR Markers Reveal Novel Potential Pathways of Chemical Stress Responses

The present de novo RNA-Seq approach yielded a *L. perenne* xenobiotic-dependent genomic database (Figures 1, 2) that was coherent with recent transcriptomic analyses of *L. perenne* (Farrell et al., 2014; Duhoux et al., 2015) and that was very useful to explore relationships between metabolic changes, stress adjustments, regulation of metabolism genes and signaling pathways under conditions of xenobiotic stress. The experimental conditions of short-term exposure and low levels of xenobiotics led, in all of the cases (glyphosate, tebuconazole, GT mixture), to no observable adverse effects, especially on root growth, which is particularly xenobiotic-sensitive (Serra et al., 2013, 2015). However, these conditions were sufficient to induce significant molecular and metabolic changes (Figures 3–6), which could be expected to reflect damage, and probably primary events in the xenobiotic response. The strong impact of these various xenobiotic treatments on the proportions of transcript functional categories and on individual gene expression established that low-level xenobiotic exposure strongly interacted with molecular regulations of gene expression. Moreover, differential effects of the different xenobiotic treatments showed that there were some xenobiotic-dependent specificities in the potential molecular mechanisms of xenobiotic action and of plant responses to xenobiotics, and also that there were specific responses for mixture (Figures 3–5). However, this global analysis also revealed the involvement of common mechanisms that are interdependent with each other and closely connected. Finally, detection of major metabolic and molecular rearrangements in leaves, which were not the primary site of exposure, demonstrated that root exposure
<figure><img src="image_10.png" /><figcaption>**FIGURE 5** | Effects of xenobiotics on metabolite levels in *L. perenne* leaves. *L. perenne* plants were grown for 7 days in the absence of pesticide and then transferred to fresh growth medium containing the different xenobiotic treatments [glyphosate (G, 1 µM), tebuconazole (T, 4 µM) and a combination of glyphosate and tebuconazole (GT, 1 µM and 4 µM, respectively)]. Metabolite levels (mean value ± SEM) were quantified after 4 days of treatment. Stars indicated statistical significance of differences between treatment and Control (p < 0.05) carried out using the Mann-Whitney-Wilcoxon test.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**FIGURE 6** | Metabolic integration of low-intensity xenobiotic stress in *L. perenne*. Black names represent quantified metabolite and gray names represent non-quantified metabolites. 2PG, 2-Phosphoglycerate; 3PG, 3-Phosphoglycerate; APS, adenosine 5'-phosphosulfate; G3P, glyceraldehyde-3-phosphate; PEP, Phosphoenolpyruvate; SAMDC, S-adenosylmethionine decarboxylase; Spd, spermidine; Spm, spermine. The red or green colors respectively indicate increase or decrease of metabolite levels in the presence of glyphosate (G), tebuconazole (T) and a combination of glyphosate and tebuconazole (GT) compared to the control condition. Solid line indicate direct relation between metabolites, and dotted line show indirect reaction.</figcaption></figure>

involved strong xenobiotic-related root-shoot communication, through xenobiotic, metabolite or signal transport, underlining the impact of such exposure.

Major genes that were found to be DE under xenobiotic exposure did not belong to classical xenobiotic or herbicide stress response pathways such as detoxification pathways involving glutathione S-transferase, cytochrome P450 or glycosyltransferase enzyme activities (Duhoux and Délye, 2013; Délye, 2013), thus suggesting that novel molecular mechanisms leading to xenobiotic-induced metabolic rearrangements have been identified. Most of the 69 DE genes, underlined without a *priori* by the DESeq analysis, surprisingly formed a coherent set as shown by the several associated pathways described below. A significant number of xenobiotic-regulated genes were linked to signal metabolisms and to signal transduction pathways, thus underlining the importance of homeostatic mechanisms and crosstalks between metabolites, carbohydrates, phytohormones, plastid-to-nucleus and light signaling under conditions of xenobiotic stress. The DESeq analysis also suggested that activities of phosphorylation linked to signal transduction

could be affected by differential expression of genes related to MAPK cascade and to ATP synthesis or hydrolysis (Table 3). Increase of phenylalanine levels in leaves in response to GT and T treatments was associated with repression of genes related to phenylalanine ammonia lyase (PAL) activity (Figures 4–6). PALs catalyze the first step of the phenylpropanoid pathway and have overlapping roles in plant growth, development, and responses to environmental stresses (Huang et al., 2010). When added at high and toxic levels, the herbicides paraquat and glyphosate induce PAL activity in plants, thus suggesting involvement of PAL in regulation of the phenylpropanoid pathway in response to oxidative stress produced by lethal herbicide levels (Duke et al., 1980; Lee et al., 2003). In *L. perenne*, low-intensity xenobiotic stress also had significant gene expression effects (Table 3, Figure 4) on the gene annotated as Cytochrome P450 reductase protein, which functions in electron transfer from NADPH to cytochrome P450 oxidases, and on a caffeate O-methyltransferase-homologous gene. These enzyme activities may be both involved in the general phenylpropanoid pathway, as described by Sundin et al. (2014) and Do et al.
(2007) for Arabidopsis, strengthening the involvement of PAL in chemical stress responses. PAL, ACC synthase, which catalyzes the synthesis of ACC for ethylene production, and auxin biosynthesis are thought to be co-regulated (Duke et al., 1980; Soeno et al., 2010). Tebuconazole treatment increased the expression of the gene annotated as β-1,3-endoglucanase and of an ACC-oxidase-annotated gene (Table 3), which functions are both related to ethylene action through regulation or metabolism (Zhong and Burns, 2003; Soeno et al., 2010; Jafari et al., 2013). These molecular markers confirmed the relationships between xenobiotic effects, phytohormone, and more particularly ethylene, dynamics and signaling, as previously reported in studies coupling hormone-signaling mutants and xenobiotic-induced stress (Sulmon et al., 2007; Weisman et al., 2010).

GO term enrichments allowed to detect by reads counting the overexpression of unigenes whose functions were related to "magnesium chelatase activity" (Supplemental Table 2). Nevertheless, comp7278_c0_seq1, whose direct function, according to its annotation, is a magnesium chelatase activity (Magnesium-chelatase subunit ChlH), presented an expression significantly affected by glyphosate and tebuconazole (Table 3, Figure 4), suggesting that overexpressed unigenes may be involved in the negative regulation of expression of the magnesium chelatase activity. The Fluorescent In Blue Light (FLU)-annotated gene was also significantly affected by glyphosate and tebuconazole single treatments (Table 3, Figure 4). FLU and ChlH are involved in regulation of tetrapyrrole biosynthesis (Kauss et al., 2012; Apitz et al., 2014) and in one of the multiple signaling pathways which mediate plastid signals controlling expression of Photosynthesis-Associated Nuclear Genes (PhANG; Mochizuki et al., 2008; Moulin et al., 2008; Zhang et al., 2013). Moreover, tetrapyrrole gene expression and ABA signaling are interconnected in chloroplast-to-nucleus signaling pathways (Voigt et al., 2010). Zhang et al. (2011) showed that resistance to concomitant high-light and herbicide stress requires the activity of ABI4, an ABA-regulated Apetala 2 (AP2)-type transcription factor, which is also involved in the coordination of several pathways such as sugar, redox, and hormonal [ABA, jasmonate (JA) and salicylic acid] signaling during stress response (Foyer et al., 2012). ABA catabolism markers were also reported to be affected in response to chemical stress treatments as glyphosate (Serra et al., 2013) and triazole including tebuconazole (Saito et al., 2006; Serra et al., 2013). Moreover, triazoles have also been shown to be inhibitors of brassinosteroid synthesis (Kaschani and van der Hoorn, 2007), of the plant sterol pathway in Sorghum (Lamb et al., 2001) and of gibberellin accumulation in rapeseed (Child et al., 1993). Thus, different hormone pathways and signaling seem to be important in numerous xenobiotic responses (Couée et al., 2013). As an example, promoter-level induction by glyphosate and hormones has been described in rice for a glutathione S-transferase involved in xenobiotic detoxification (Hu et al., 2011). However, further experimental studies need to be carried out in order to clearly establish the role of these hormones in plant chemical stress responses.

## Low Levels of Xenobiotics Have Significant Effects on Multiple Low Energy and Low Carbohydrate Regulations

Genes related to components of photosynthetic reaction centers presented patterns of expression that were largely influenced by xenobiotics. The gene annotated as photosystem II reaction center protein D presented the highest observed induction by G and T and repression by GT. More generally, differential regulation of genes related to photosynthesis in the presence of sub-lethal concentrations of xenobiotics was specific of the xenobiotic involved (Table 3). Some herbicides, such as atrazine, act through binding to D1 protein of photosystem II (PSII) reaction center, thus blocking electron transfer to the plastoquinone pool (Rutherford and Krieger-Liszkay, 2001; Sulmon et al., 2004). However, neither glyphosate nor tebuconazole are known to affect photosynthesis through this mechanism of action. These variations may thus originate from non-target effects, in accordance with the maintenance of aromatic amino acid levels under G treatment (Figure 5), suggesting that the canonical EPSPS target of glyphosate was not affected. These differential regulations of photosynthesis-related genes by chemical treatments did not however result in modifications of photosynthetic parameters (Serra et al., 2015), as also reported by Ramel et al. (2007) and Faus et al. (2015), underlining adjustment processes to maintain, at least temporarily, efficient carbon assimilation.

Variations of carbohydrate-related and energy-related GO term categories (Figure 3) and, more directly, the variations of soluble sugars (Figure 5) strongly indicated that carbohydrate and energy metabolisms were involved in the responses to low-intensity xenobiotic stress. The level of water-soluble carbohydrates is particularly important as a primary source of energy in *L. perenne* (Smith et al., 1998). Some of the metabolic changes induced by xenobiotics, especially soluble sugar depletion and lactate accumulation (Figures 5, 6), therefore suggested a situation of energy imbalance. All of this strongly suggested that low-intensity xenobiotic stress could be perceived as a central energy perturbation. It must be highlighted that such a relationship has been described in other plant-xenobiotic experimental systems (Qian et al., 2014; Serra et al., 2015). Stress-induced energy deficit typically leads to the activation of genes involved in catabolism (proteolysis, amino acid catabolism, sugar degradation, lipid mobilization) and repression of genes involved in protein synthesis (Valluru and Van den Ende, 2011; Ramel et al., 2012; Serra et al., 2013). The strong variations of amino acid levels (Figure 5) and of amino acid-related genes (Figure 4, Table 3) may reflect the mobilization of amino acid substrates (Valluru and Van den Ende, 2011; Ramel et al., 2012; Serra et al., 2013; Duhoux et al., 2015). Further studies are needed to characterize the relations between amino acid and carbohydrate metabolisms, energy perturbation and this particular NOAE situation produced by short exposure to low doses of xenobiotics.

The potential variations of spermine and spermidine levels, that could be deduced from the strong increase in unigenes related to the "spermine biosynthetic process" category (Figure 3A), and from the decrease of putrescine levels in
shoots from xenobiotic-treated plants (Figure 5), also suggested that low-level xenobiotic exposure caused rearrangements of N-related metabolisms.

Situations of low carbohydrate and low energy can be strongly associated with oxidative stress processes (Couée et al., 2006; Baena-González et al., 2007; Dietrich et al., 2011; Valluru and Van den Ende, 2011). The differential expression of genes linked to the pentose phosphate shunt (Table 3) could thus be related to variations of carbohydrate levels and to the supply of reducing equivalents for antioxidant defenses (Couée et al., 2006). Low-intensity xenobiotic stress modified the expression of other genes related to ROS defense. A Glycolate oxidase (GOX1)-annotated gene was found to be regulated in a contrasted manner by glyphosate, tebuconazole and GT treatments (Table 3). Glycolate oxidase is involved in the oxidation of glycolate leading to the production of glyoxylate and H₂O₂ during photorespiration (Fahnenstich et al., 2008). In parallel, the levels of photorespiration metabolites (Figures 5, 6) and the expression of a Glycine Decarboxylase-annotated gene (Figure 4) were strongly affected by xenobiotic treatments. These results suggested that these conditions of mild chemical stress, and in relation with carbohydrate, energy and ROS processes, could affect photorespiration.

## Various Molecular Markers Point Out to the Importance of SnRKs in Xenobiotic Stress Responses

Sugar and energy limitation induced by abiotic stresses are regulated in plants by an array of metabolic sensors (Baena-González and Sheen, 2008; Dietrich et al., 2011). One of these sensors, SnRK1, is involved in responses to a range of stresses that limit photosynthesis and respiration, including the PSII-inhibiting herbicide, DCMU (Baena-González et al., 2007). Sugars, and more particularly Glc-6-P produced by hexokinase (HXK) and Tre-6-P, can block SnRK1 activity (Baena-González and Sheen, 2008; Ramon et al., 2008; Valluru and Van den Ende, 2011; Nunes et al., 2013). Our previous global analysis in *L. perenne*, combining leaves and root metabolites after short exposure and direct and long exposure, showed that various xenobiotic treatments mostly increased Glc-6-P and its correlated co-metabolite, Fru-6-P, levels (Serra et al., 2015). Under the present conditions of short exposure, Glc-6-P remained relatively stable in leaves, suggesting that SnRK1 may not be inhibited by the phosphorylated sugar.

In the present case of *L. perenne*-xenobiotic interactions, several genes, annotated as F-type H⁺-transporting ATPase, H⁺-ATPase, and ATP synthase, were DE in the presence of xenobiotics (Table 3). Moreover, repression of a gene potentially involved in ATP hydrolysis and induction of ATP synthesis genes suggested some kind of cellular mobilization for energy restoration under conditions of xenobiotic exposure (Table 3). Decrease of, or threat to, ATP levels occurring under environmental stresses, including xenobiotic stress, were also related to the involvement of mammalian homolog of SnRK1, the AMP-activated kinase (AMPK) energy sensor (Blättler et al., 2007; Hardie et al., 2012).

SnRK1 and several other serine/threonine kinases, especially SnRK2 and SnRK3, have activities that link stress and ABA signaling with metabolic signaling (Halford and Hey, 2009; Yu et al., 2014). In Arabidopsis, 38 members of the SnRK family have been described. The different SnRK subfamilies, SnRK1, SnRK2 and SnRK3, have respectively 3, 10 and 25 representatives (Halford and Hey, 2009). All of these SnRKs are associated with metabolic regulations and stress responses. The present transcriptome reveals that *L. perenne* possesses 146 unigenes with high homologies (e < 10e-6) with these SnRK subfamilies. In comparison, a recent study has identified 44 SnRKs in *Brachypodium distachyon* that respond to multiple stresses and stress-related signal molecules like ABA, ethylene and H₂O₂ (Wang et al., 2015). At least some of *L. perenne* SnRK genes were strongly modulated by low-intensity xenobiotic stress (Table 3, Figure 4, Supplemental Figure 1), thus strengthening the above-discussed link with the involvement of ABA signaling in xenobiotic responses. ABA regulates the activity of SnRK2 and SnRK3 subfamilies which are part of signaling pathways involved in responses to abiotic stresses such as nutrient limitation, drought, cold, salt, and osmotic stresses (Coello et al., 2011). It has also been shown that SnRK2 could be activated by Ca²⁺ (Coello et al., 2012). Among SnRK2-interacting proteins, Protein Phosphatase 2Cs (PP2Cs) act as negative modulators in relation with ABA signaling (Fujita et al., 2013; Danquah et al., 2014). Interestingly, a PP2C-annotated gene was found to be regulated by xenobiotic treatments, at least in the DESeq analysis (Table 3). Moreover, the SnRK2.3-annotated unigene was DE after xenobiotic treatments (Figure 4). DESeq and qRT-PCR analysis also showed that a transcript related to SnRK3.12 (CIPK9) was repressed by G treatment, which, in contrast, induced the expression of a calcium sensing receptor gene and of an aldehyde dehydrogenase gene. Such aldehyde dehydrogenases have been involved in calcium-mediated signaling and stress-regulated detoxification in young plants (Stiti et al., 2011). SnRK3s interact with calcineurin B-like (CBL) calcium-binding proteins, hence their description as CBL-interacting kinases (CIPKs; Coello et al., 2011), and these CIPKs carry out crosstalks between Ca²⁺ and ABA signaling in the responses to abiotic stresses (Yu et al., 2014). Despite some contradictory studies about the sensitivity of plants lacking *CIPK9* on low potassium media, a function in potassium homeostasis has also been suggested (Pandey et al., 2007; Liu et al., 2013b; Yu et al., 2014). All of these data therefore suggest that SnRKs and interacting proteins could be implicated in xenobiotic responses, cell homeostasis under xenobiotic stress and xenobiotic stress adaptation through complex regulations that involve hormone signaling, calcium signaling and potassium homeostasis. Further studies, using hormone inhibitors or antagonists and mutants related to SnRK or hormone signaling, would allow to understand the role of SnRK-related signaling in plant responses to xenobiotics. Interestingly, herbicide resistance in the agricultural weed *Ipomoea purpurea* has been associated with the differential expression of genes involved in such cellular signaling, with high glyphosate treatment inducing an AtPK7-like serine/threonine-protein kinase and repressing genes encoding receptor-like kinases (Leslie and Baucom, 2014).
CONCLUSION

The present results are in line with the involvement of genes
linked to SnRKs and other protein kinases and to hormonal
regulation in response to low-intensity chemical stresses in L.
*perenne*. Our results emphasize that responses to xenobiotic
stresses require signal transduction mechanisms involving
protein-kinases and protein-phosphatases, in accordance with
other studies (Sulmon et al., 2007; Fukudome et al., 2014). Our
parallel metabolic and molecular approaches establish a vision of
responses to mild chemical stress potentially involving complex
and interconnected regulation processes. Energy dysfunction is
likely to be related to carbohydrate dynamics and regulation.
The present results therefore suggested the existence of complex
relations between sugar signaling, the signaling pathways
associated with phytohormones and calcium (León and Sheen,
2003; Rolland et al., 2006; Kudla et al., 2010; Biswal et al., 2011)
and the molecular responses to organic xenobiotics (Ramel et al.,
2012; Serra et al., 2013). Finally the close relations between
protein-kinase regulations and ABA signaling suggest hormonal
actions in response to mild chemical stress. The exact role of
these processes in xenobiotic responses remains to be elucidated.
This role is likely to be connected to molecular and physiological
processes involved in plant growth and development. However,
the nature of the primary events leading to such signaling is
still completely unknown. Interestingly, herbicide resistance in
the agricultural weed *Ipomoea purpurea* has been associated
with the differential expression of genes involved in such
cellular signaling, with high glyphosate treatment inducing
an AtPK7-like serine/threonine-protein kinase and repressing
genes encoding receptor-like kinases (Leslie and Baucom, 2014).
Situations of NTSR observed in rye-grasses have been correlated
with mechanisms of gene expression regulation (Duhoux et al.,

2015; Salas et al., 2015). Identification and characterization of key
steps of the corresponding signal transduction processes should
therefore be of particular interest for herbicide management in
crop systems and for restoration-phytoremediation activities in
ecological engineering.

ACKNOWLEDGMENTS

We are grateful to Dr Erwan Corre from the CNRS-UPMC
ABIMS bioinformatics platform (Station biologique de Roscoff,
France) for providing help and support with bioinformatics
analysis, the ABGC technical core facility (UMR CNRS 6553,
University of Rennes 1, France) for biochemical analyses and the
Biogenouest Environmental and Functional Genomics Platform
core facility (UMS CNRS 3343, Observatoire des Sciences de
l'Univers de Rennes, France) for transcriptome sequencing.
We also wish to thank Jean-Luc Foulon (UMR CNRS 6553,
University of Rennes 1, France) for help with plant culture
systems and Sandra Rigaud (UMR CNRS 6553, University of
Rennes 1, France) for administrative and accounting work.
This work was supported by the interdisciplinary program
"Ingénierie écologique" from the Centre National de la Recherche
Scientifique (CNRS, France) and by the Fondation pour la
Recherche sur la Biodiversité (FRB, France). AS is supported
by a doctoral scholarship from the Brittany regional council
(France).

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found
online at: http://journal.frontiersin.org/article/10.3389/fpls.2015.
01124

REFERENCES

Aliferis, K. A., Materzok, S., Paziotou, G. N., and Chrysayi-Tokousbalides, M. (2009). *Lemna minor* L. as a model organism for ecotoxicological studies performing ¹H NMR fingerprinting. *Chemosphere* **76**, 967–973. doi: **10.1016/j.chemosphere.2009.04.025**

Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., Miller, W., et al. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. *Nucleic Acids Res.* **25**, 3389–3402. doi: 10.1093/nar/25.1.7.3389

Anders, S., McCarthy, D. J., Chen, Y., Okoniewski, M., Smyth, G. K., Huber, W., et al. (2013). Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. *Nat. Protoc.* **8**, 1765–1786. doi: **10.1038/nprot.2013.099**

Anjum, N. A., Gill, R., Kaushik, M., Hasanuzzaman, M., Pereira, E., Ahmad, I., et al. (2015). ATP-sulfurylase, sulfur-compounds, and plant stress tolerance. *Front. Plant Sci.* **6**:210. doi: 10.3389/fpls.2015.00210

Apitz, J., Schmied, J., Lehmann, M. J., Hedtke, B., and Grimm, B. (2014). GluTR2 complements a *hema1* mutant lacking glutamyl-tRNA reductase 1, but is differently regulated at the post-translational level. *Plant Cell Physiol.* **55**, 645–657. doi: 10.1093/pcp/pcu016

Arias-Estévez, M., López-Periago, E., Martínez-Carballo, E., Simal-Gándara, J., Mejuto, J.-C., and García-Río, L. (2008). The mobility and degradation of pesticides in soils and the pollution of groundwater resources. *Agric. Ecosyst. Environ.* **123**, 247–260. doi: 10.1016/j.agee.2007.07.011

Baena-González, E., Rolland, F., Thevelein, J. M., and Sheen, J. (2007). A central integrator of transcription networks in plant stress and energy signalling. *Nature* **448**, 938–943. doi: 10.1038/nature06069

Baena-González, E., and Sheen, J. (2008). Convergent energy and stress signaling. *Trends Plant Sci.* **13**, 474–482. doi: 10.1016/j.tplants.2008.06.006

Barbehenn, R. V., Karowe, D. N., and Spickard, A. (2004). Effects of elevated atmospheric CO₂ on the nutritional ecology of C3 and C4 grass-feeding caterpillars. *Oecologia* **140**, 86–95. doi: 10.1007/s00442-004-1572-9

Belz, R. G., and Duke, S. O. (2014). Herbicides and plant hormesis. *Pest Manag. Sci.* **70**, 698–707. doi: 10.1002/ps.3726

Bidar, G., Pruvot, C., Garçon, G., Verdin, A., Shirali, P., and Douay, F. (2009). Seasonal and annual variations of metal uptake, bioaccumulation, and toxicity in *Trifolium repens* and *Lolium perenne* growing in a heavy metal-contaminated field. *Environ. Sci. Pollut. Res. Int.* **16**, 42–53. doi: 10.1007/s11356-008-0021-4

Biswal, B., Joshi, P. N., Raval, M. K., and Biswal, U. C. (2011). Photosynthesis, a global sensor of environmental stress in green plants: stress signalling and adaptation. *Curr. Sci.* **101**, 47–56.

Blättler, S. M., Rencurel, F., Kaufmann, M. R., and Meyer, U. A. (2007). In the regulation of cytochrome P450 genes, phenobarbital targets LKB1 for necessary activation of AMP-activated protein kinase. *Proc. Natl. Acad. Sci. U.S.A.* **104**, 1045–1050. doi: 10.1073/pnas.0610216104

Casler, M. D., and Duncan, R. R. (2003). *Turfgrass Biology, Genetics, and Breeding*. Hoboken, NJ: John Wiley & Sons, Inc.

Cedergreen, N. (2014). Quantifying synergy: a systematic review of mixture toxicity studies within environmental toxicology. *PLoS ONE* **9:e96580**. doi: **10.1371/journal.pone.0096580**
Child, R. D., Evans, D. E., Allen, J., and Arnold, G. M. (1993). Growth responses in oilseed rape (*Brassica napus* L.) to combined applications of the triazole chemicals triapenthenol and tebuconazole and interactions with gibberellin. *Plant Growth Regul.* 13, 203-212. doi: 10.1007/BF00024263

Coello, P., Hey, S. J., and Halford, N. G. (2011). The sucrose non-fermenting-1-related (SnRK) family of protein kinases: potential for manipulation to improve stress tolerance and increase yield. *J. Exp. Bot.* 62, 883-893. doi: 10.1093/jxb/erq331

Coello, P., Hirano, E., Hey, S. J., Muttucumaru, N., Martinez-Barajas, E., Parry, M. A., et al. (2012). Evidence that abscisic acid promotes degradation of SNF1-related protein kinase (SnRK) 1 in wheat and activation of a putative calcium-dependent SnRK2. *J. Exp. Bot.* 63, 913-924. doi: 10.1093/jxb/err320

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. *Bioinformatics* 21, 3674-3676. doi: 10.1093/bioinformatics/bti610

Costantini, D., Metcalfe, N. B., and Monaghan, P. (2010). Ecological processes in a hermetic framework. *Ecol. Lett.* 13, 1435-1447. doi: 10.1111/j.1461-0248.2010.01531.x

Couée, I., Serra, A. A., Ramel, F., Gouesbet, G., and Sulmon, C. (2013). Physiology and toxicology of hormone-disrupting chemicals in higher plants. *Plant Cell Rep.* 32, 933-941. doi: 10.1007/s00299-013-1428-z

Couée, I., Sulmon, C., Gouesbet, G., and El Amrani, A. (2006). Involvement of soluble sugars in reactive oxygen species balance and responses to oxidative stress in plants. *J. Exp. Bot.* 57, 449-459. doi: 10.1093/jxb/erj027

D'Orazio, V., Ghanem, A., and Senesi, N. (2013). Phytoremediation of pyrene contaminated soils by different plant species. *CLEAN* 41, 377-382. doi: 10.1002/clen.201100653

Danquah, A., de Zelicourt, A., Colcombet, J., and Hirt, H. (2014). The role of ABA and MAPK signaling pathways in plant abiotic stress responses. *Biotechnol. Adv.* 32, 40-52. doi: 10.1016/j.biotechadv.2013.09.006

Das, M., Reichman, J. R., Haberer, G., Welzl, G., Aceituno, F. F., Mader, M. T., et al. (2010). A composite transcriptional signature differentiates responses towards closely related herbicides in *Arabidopsis thaliana* and *Brassica napus*. *Plant Mol. Biol.* 72, 545-556. doi: 10.1007/s11103-009-9590-y

Dear, B. S., Sandral, G. A., and Wilson, B. C. D. (2006). Tolerance of perennial pasture grass seedlings to pre- and post-emergent grass herbicides. *Aust. J. Exp. Agr.* 46, 637-644. doi: 10.1071/EA04173

Délye, C. (2013). Unravelling the genetic bases of non-target-site-based resistance (NTSR) to herbicides: a major challenge for weed science in the forthcoming decade. *Pest Manag. Sci.* 69, 176-187. doi: 10.1002/ps.3318

Dévérier, M. H., Mazellier, P., Aït-Aïssa, S., and Budzinski, H. (2011). New challenges in environmental analytical chemistry: identification of toxic compounds in complex mixtures. *Comptes Rendus Chimie* 14, 766-779. doi: 10.1016/j.crci.2011.04.006

Diekmann, K., Hodkinson, T. R., Wolfe, K. H., van den Bekerom, R., Dix, P. J., and Barth, S. (2009). Complete chloroplast genome sequence of a major allogamous forage species, Perennial Ryegrass (*Lolium perenne* L.). *DNA Res.* 16, 165-176. doi: 10.1093/dnares/dsp008

Dietrich, K., Weltmeier, F., Ehlert, A., Weiste, C., Stahl, M., Harter, K., et al. (2011). Heterodimers of the Arabidopsis transcription factors bZIP1 and bZIP53 reprogram amino acid metabolism during low energy stress. *Plant Cell* 23, 381-395. doi: 10.1105/tpc.110.075390

Do, C. T., Pollet, B., Thévenin, J., Sibout, R., Denoue, D., Barrière, Y., et al. (2007). Both caffeoyl Coenzyme A 3-O-methyltransferase 1 and caffeic acid O-methyltransferase 1 are involved in redundant functions for lignin, flavonoids and sinapoyl malate biosynthesis in Arabidopsis. *Planta* 226, 1117-1129. doi: 10.1007/s00425-007-0558-3

Duhoux, A., Carrère, S., Gouzy, J., Bonin, L., and Délye, C. (2015). RNA-Seq analysis of rye-grass transcriptomic response to an herbicide inhibiting acetolactate-synthase identifies transcripts linked to non-target-site-based resistance. *Plant Mol. Biol.* 87, 473-487. doi: 10.1007/s11103-015-0292-3

Duhoux, A., and Délye, C. (2013). Reference genes to study herbicide stress response in *Lolium* sp.: up-regulation of P450 genes in plants resistant to acetolactate-synthase inhibitors. *PLoS ONE* 8:e63576. doi: 10.1371/journal.pone.0063576

Duke, S. O., Hoagland, R. E., and Elmore, C. D. (1980). Effects of glyphosate on metabolism of phenolic compounds: V. l-alpha-aminoxy-beta-phenylpropionic acid and glyphosate effects on phenylalanine ammonia-lyase in soybean seedlings. *Plant Physiol.* 65, 17-21. doi: 10.1104/pp.65.1.17

Duke, S. O., and Powles, S. B. (2008). Glyphosate: a once-in-a-century herbicide. *Pest Manag. Sci.* 64, 319-325. doi: 10.1002/ps.1518

Fahnenstich, H., Scarpeci, T. E., Valle, E. M., Flügge, U. I., and Maurino, V. G. (2008). Generation of hydrogen peroxide in chloroplasts of Arabidopsis overexpressing glycolate oxidase as an inducible system to study oxidative stress. *Plant Physiol.* 148, 719-729. doi: 10.1104/pp.108.126789

Farrell, J. D., Byrne, S., Paina, C., and Asp, T. (2014). *De novo* assembly of the perennial ryegrass transcriptome using an RNA-seq strategy. *PLoS ONE* 9:e103567. doi: 10.1371/journal.pone.0103567

Faus, I., Zabalza, A., Santiago, J., Nebauer, S. G., Royuela, M., Serrano, R., et al. (2015). Protein kinase GCN2 mediates responses to glyphosate in Arabidopsis. *BMC Plant Biol.* 15:14-25. doi: 10.1186/s12870-014-0378-0

Fiil, A., Lenk, I., Petersen, K., Jensen, C. S., Nielsen, K. K., Schejbel, B., et al. (2011). Nucleotide diversity and linkage disequilibrium of nine genes with putative effects on flowering time in perennial ryegrass (*Lolium perenne* L.). *Plant Sci.* 180, 228-237. doi: 10.1016/j.plantsci.2010.08.015

Foyer, C. H., Kerchev, P. I., and Hancock, R. D. (2012). The ABA-INSENSITIVE-4 (ABI4) transcription factor links redox, hormone and sugar signaling pathways. *Plant Signal Behav.* 7, 276-281. doi: 10.4161/psb.18770

Frankart, C., Eullaffroy, P., and Vernet, G. (2002). Photosynthetic responses of *Lemna minor* exposed to xenobiotics, copper, and their combinations. *Ecotoxico. Environ. Saf.* 53, 439-445. doi: 10.1016/S0147-6513(02)00003-9

Frye, C. A., Bo, E., Calamandrei, G., Calzà, L., Dessi-Fulgheri, F., Fernández, M., et al. (2011). Endocrine disrupters: a review of some sources, effects, and mechanisms of actions on behaviour and neuroendocrine systems. *J. Neuroendocrinol.* 24, 144-159. doi: 10.1111/j.1365-2826.2011.02229.x

Fujita, Y., Yoshida, T., and Yamaguchi-Shinozaki, K. (2013). Pivotal role of the AREB/ABF-SnRK2 pathway in ABRE-mediated transcription in response to osmotic stress in plants. *Physiol Plant.* 147, 15-27. doi: 10.1111/j.1399-3054.2012.01635.x

Fukudome, A., Aksoy, E., Wu, X., Kumar, K., Jeong, I. S., May, K., et al. (2014). Arabidopsis CPL4 is an essential C-terminal domain phosphatase that suppresses xenobiotic stress responses. *Plant J.* 80, 27-39. doi: 10.1111/tpj.12612

Furet, P. M., Berthier, A., Decau, M. L., Morvan-Bertrand, A., Prud'homme, M. P., Noiraud-Romy, N., et al. (2012). Differential regulation of two sucrose transporters by defoliation and light conditions in perennial ryegrass. *Plant Physiol. Biochem.* 61, 88-96. doi: 10.1016/j.plaphy.2012.09.011

Galston, A. W., and Sawhney, R. K. (1990). Polyamines in plant physiology. *Plant Physiol.* 94:406-410. doi: 10.1104/pp.94.2.406

Goecks, J., Nekrutenko, A., Taylor, J., and The Galaxy Team (2010). Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. *Genome Biol.* 11:R86. doi: 10.1186/gb-2010-11-8-r86

Gomes, M. P., Smedbol, E., Chalifour, A., Hénault-Ethier, L., Labrecque, M., Lepage, L., et al. (2014). Alteration of plant physiology by glyphosate and its by-product aminomethylphosphonic acid: an overview. *J. Exp. Bot.* 65, 4691-4703. doi: 10.1093/jxb/eru269

Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2013). *De novo* transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. *Nat. Protoc.* 8, 1494-1512. doi: 10.1038/nprot.2013.084

Halford, N. G., and Hey, S. J. (2009). SNF1-related protein kinases (SnRKs) act within an intricate network that links metabolic and stress signalling in plants. *Biochem. J.* 419, 247-259. doi: 10.1042/BJ20082408

Hardie, D. G., Ross, F. A., and Hawley, S. A. (2012). AMPK: a nutrient and energy sensor that maintains energy homeostasis. *Nat. Rev. Mol. Cell Biol.* 13, 251-262. doi: 10.1038/nrm3111

Helander, M., Saloniemi, I., and Saikkonen, K. (2012). Glyphosate in northern ecosystems. *Trends Plant Sci.* 17, 569-574. doi: 10.1016/j.tplants.2012.05.008

Hernández, A. F., Parrón, T., Tsatsakis, A. M., Requena, M., Alarcón, R., and López-Guarnido, O. (2012). Toxic effects of pesticide mixtures at a molecular level: their relevance to human health. *Toxicology* 307, 136-145. doi: 10.1016/j.tox.2012.06.009
Hertzberg, R. C., and MacDonell, M. M. (2002). Synergy and other ineffective mixture risk definitions. *Sci. Total Environ.* 288, 31–42. doi: 10.1016/S0048-9697(01)01113-5

Hu, T., He, S., Yang, G., Zeng, H., Wang, G., Chen, Z., et al. (2011). Isolation and characterization of a rice glutathione S-transferase gene promoter regulated by herbicides and hormones. *Plant Cell Rep.* 30, 539–549. doi: 10.1007/s00299-010-0964-z

Huang, H. H., Xu, L. L., Tong, Z. K., Lin, E. P., Liu, Q. P., Cheng, L. J., et al. (2012). *De novo* characterization of the Chinese fir (*Cunninghamia lanceolata*) transcriptome and analysis of candidate genes involved in cellulose and lignin biosynthesis. *BMC Genomics* 13:648. doi: 10.1186/1471-2164-13-648

Huang, J., Gu, M., Lai, Z., Fan, B., Shi, K., Zhou, Y. H., et al. (2010). Functional analysis of the Arabidopsis PAL gene family in plant growth, development, and response to environmental stress. *Plant Physiol.* 153, 1526–1538. doi: 10.1104/pp.110.157370

Islam, M. S., Studer, B., Byrne, S. L., Farrell, J. D., Panitz, F., Bendixen, C., et al. (2013). The genome and transcriptome of perennial ryegrass mitochondria. *BMC Genomics* 14:202. doi: 10.1186/1471-2164-14-202

Ivanov, S., Shopova, E., Kerchev, P., Sergiev, I., Miteva, L., Polizov, D., et al. (2013). Long-term impact of sublethal atrazine perturbs the redox homeostasis in pea (*Pisum sativum* L.) plants. *Protoplasma* 250, 95–102. doi: 10.1007/s00709-012-0378-6

Jafari, Z., Haddad, R., Hosseini, R., and Garoosi, G. (2013). Cloning, identification and expression analysis of ACC oxidase gene involved in ethylene production pathway. *Mol. Biol. Rep.* 40, 1341–1350. doi: 10.1007/s11033-012-2178-7

Kaschani, F., and van der Hoorn, R. (2007). Small molecule approaches in plants. *Curr. Opin. Chem. Biol.* 11, 88–98. doi: 10.1016/j.cbpa.2006.11.038

Kauss, D., Bischof, S., Steiner, S., Apel, K., and Meskauskiene, R. (2012). FLU, a negative feedback regulator of tetrapyrrole biosynthesis, is physically linked to the final steps of the Mg<sup>++</sup>-branch of this pathway. *FEBS Lett.* 586, 211–216. doi: 10.1016/j.febset.2011.12.029

Kertesz, N., Samson, J., Debacker, C., Wu, H., and Labastie, M. C. (2002). Cloning and characterization of human and mouse SNRK sucrose non-fermenting protein (SNF-1)-related kinases. *Gene* 294, 13–24. doi: 10.1016/S0378-1119(02)00829-6

Köhler, H. R., and Triebskorn, R. (2013). Wildlife ecotoxicology of pesticides: can we track effects to the population level and beyond? *Science* 341, 759–765. doi: 10.1126/science.1237591

Kranner, I., Minibayeva, F. V., Beckett, R. P., and Seal, C. E. (2010). What is stress? Concepts, definitions and applications in seed science. *New Phytol.* 188, 655–673. doi: 10.1111/j.1469-8137.2010.03461.x

Kudla, J., Batistic, O., and Hashimoto, K. (2010). Calcium signals: the lead currency of plant information processing. *Plant Cell* 22, 541–563. doi: 10.1105/tpc.109.072686

Lamb, D. C., Cannieux, M., Warrilow, A. G., Bak, S., Kahn, R. A., Manning, N. J., et al. (2001). Plant sterol 14α-demethylase affinity for azole fungicides. *Biochem. Biophys. Res. Commun.* 284, 845–849. doi: 10.1006/bbrc.2001.5010

Lee, B. K., Park, M. R., Srinivas, B., Chun, J. C., Kwon, I. S., Chung, I. M., et al. (2003). Induction of phenylalanine ammonia-lyase gene expression by paraquat and stress-related hormones in *Rehmannia glutinosa*. *Mol. Cells* 16, 34–39.

León, P., and Sheen, J. (2003). Sugar and hormone connections. *Trends Plant Sci.* 8, 110–116. doi: 10.1016/S1360-1385(03)00011-6

Leslie, T., and Baucom, R. S. (2014). *De novo* assembly and annotation of the transcriptome of the agricultural weed *Ipomoea purpurea* uncovers gene expression changes associated with herbicide resistance. *G3 (Bethesda)*. 4, 2035–2047. doi: 10.1534/g3.114.013508

Li, Z., Cheng, Y., Cui, J., Zhang, P., Zhao, H., and Hu, S. (2015). Comparative transcriptome analysis reveals carbohydrate and lipid metabolism blocks in *Brassica napus* L. male sterility induced by the chemical hybridization agent monosulfuron ester sodium. *BMC Genomics* 16:206. doi: 10.1186/s12864-015-1388-5

Liu, L. L., Ren, H. M., Chen, L. Q., Wang, Y., and Wu, W. H. (2013b). A protein kinase, calcineurin B-like protein-interacting protein Kinase9, interacts with calcium sensor calcineurin B-like Protein3 and regulates potassium homeostasis under low-potassium stress in Arabidopsis. *Plant Physiol.* 161, 266–277. doi: 10.1104/pp.112.206896

Liu, M., Hulting, A. G., and Mallory-Smith, C. A. (2013a). Characterization of multiple-herbicide-resistant Italian ryegrass (*Lolium perenne* spp. *multiflorum*). *Pest Manag. Sci.* 70, 1145–1150. doi: 10.1002/ps.3665

Londo, J. P., McKinney, J., Schwartz, M., Bollman, M., Sagers, C., and Watrud, L. (2014). Sub-lethal glyphosate exposure alters flowering phenology and causes transient male-sterility in *Brassica* spp. *BMC Plant Biol.* 14:70. doi: 10.1186/1471-2229-14-70

Lydy, M., Belden, J., Wheelock, C., Hammock, B., and Denton, D. (2004). Challenges in regulating pesticide mixtures. *Ecol. Soc.* 9:1.

Medina, M. H., Correa, J. A., and Barata, C. (2007). Microevolution due to pollution: possible consequences for ecosystem responses to toxic stress. *Chemosphere* 67, 2105–2114. doi: 10.1016/j.chemosphere.2006.12.024

Mochizuki, N., Tanaka, R., Tanaka, A., Masuda, T., and Nagatani, A. (2008). The steady-state level of Mg-protoporphyrin IX is not a determinant of plastid-to-nucleus signaling in Arabidopsis. *Proc. Natl. Acad. Sci. U.S.A.* 105, 15184–15189. doi: 10.1073/pnas.0803245105

Moulin, M., McCormac, A. C., Terry, M. J., and Smith, A. G. (2008). Tetrapyrrole profiling in Arabidopsis seedlings reveals that retrograde plastid nuclear signaling is not due to Mg-protoporphyrin IX accumulation. *Proc. Natl. Acad. Sci. U.S.A.* 105, 15178–15183. doi: 10.1073/pnas.0803054105

Nadar, H. M., Clegg, M. D., and Maranville, J. W. (1975). Promotion of Sorghum callus growth by the s-triazine herbicides. *Plant Physiol.* 56, 747–751. doi: 10.1104/pp.56.6.747

Nobels, I., Spanoghe, P., Haesaert, G., Robbens, J., and Blust, R. (2011). Toxicity ranking and toxic mode of action evaluation of commonly used agricultural adjuvants on the basis of bacterial gene expression profiles. *PLoS ONE* 6:e24139. doi: 10.1371/journal.pone.0024139

Nunes, C., Primavesi, L. F., Patel, M. K., Martinez-Barajas, E., Powers, S. J., Sagar, R., et al. (2013). Inhibition of SnRK1 by metabolites: tissue-dependent effects and cooperative inhibition by glucose 1-phosphate in combination with trehalose 6-phosphate. *Plant Physiol. Biochem.* 63, 89–98. doi: 10.1016/j.plaphy.2012.11.011

Pandey, G. K., Cheong, Y. H., Kim, B. G., Grant, J. J., Li, L., and Luan, S. (2007). CIPK9: a calcium sensor-interacting protein kinase required for low-potassium tolerance in *Arabidopsis*. *Cell Res.* 17, 411–421. doi: 10.1038/cr.2007.39

Patty, L., Real, B., and Gril, J. J. (1997). The use of grassed buffer strips to remove pesticides, nitrate and soluble phosphorus compounds from runoff water. *Pestic. Sci.* 49, 243–251.

Pfleeger, T., Blakely-Smith, M., King, G., Lee, E. H., Plocher, M., and Olszyk, D. (2012). The effects of glyphosate and aminopyralid on a multi-species plant field trial. *Ecotoxicology* 21, 1771–1787. doi: 10.1007/s10646-012-0912-5

Potter, T. L., Bosch, D. D., and Strickland, T. C. (2014). Comparative assessment of herbicide and fungicide runoff risk: a case study for peanut production in the Southern Atlantic Coastal Plain (USA). *Sci. Total Environ.* 490, 1–10. doi: 10.1016/j.scitotenv.2014.04.034

Pouyat, R., Yesilonis, I., and Golubiewski, N. (2009). A comparison of soil organic carbon stocks between residential turf grass and native soil. *Urban Ecosyst.* 12, 45–62. doi: 10.1007/s11252-008-0059-6

Qian, H., Tsuji, T., Endo, T., and Sato, F. (2014). PGR5 and NDH pathways in photosynthetic cyclic electron transfer respond differently to sublethal treatment with photosystem-interfering herbicides. *J. Agric. Food. Chem.* 62, 4083–4089. doi: 10.1021/jf500143f

Ramel, F., Sulmon, C., Bogard, M., Couée, I., and Gouesbet, G. (2009). Differential dynamics of reactive oxygen species and antioxidative mechanisms during atrazine injury and sucrose-induced tolerance in *Arabidopsis thaliana* seedlings. *BMC Plant Biol.* 9:28. doi: 10.1186/1471-2229-9-28

Ramel, F., Sulmon, C., Cabello-Hurtado, F., Taconnat, L., Martin-Magniette, M. L., Renou, J. P., et al. (2007). Genome-wide interacting effects of sucrose and herbicide-mediated stress in *Arabidopsis thaliana*: novel insights into atrazine toxicity and sucrose-induced tolerance. *BMC Genomics* 8:450. doi: 10.1186/1471-2164-8-450

Ramel, F., Sulmon, C., Serra, A. A., Gouesbet, G., and Couée, I. (2012). Xenobiotic sensing and signalling in higher plants. *J. Exp. Bot.* 63, 3999–4014. doi: 10.1093/jxb/ers102

Ramon, M., Rolland, F., and Sheen, J. (2008). Sugar sensing and signaling. *Arabidopsis Book* 6:e0117. doi: 10.1199/tab.0117
Ren, X., Liu, T., Dong, J., Sun, L., Yang, J., Zhu, Y., et al. (2012). Evaluating de Bruijn graph assemblers on 454 transcriptomic data. *PLoS ONE* 7:e51188. doi: 10.1371/journal.pone.0051188

Rolland, F., Baena-Gonzalez, E., and Sheen, J. (2006). Sugar sensing and signaling in plants: conserved and novel mechanisms. *Annu. Rev. Plant Biol.* 57, 675–709. doi: 10.1146/annurev.arplant.57.032905.105441

Rozen, S., and Skaletsky, H. J. (2000). "Primer3 on the WWW for general users and for biologist programmers", in *Bioinformatics Methods and Protocols: Methods in Molecular Biology*, eds S. Misener, and S. A. Krawetz (Totowa, NJ: Humana Press Inc.), 365–386.

Rutherford, A. W., and Krieger-Liszkay, A. (2001). Herbicide-induced oxidative stress in photosystem II. *Trends Biochem. Sci.* 26, 648–653. doi: 10.1016/S0968-0004(01)01953-3

Saito, S., Okamoto, M., Shinoda, S., Kushiro, T., Koshiba, T., Kamiya, Y., et al. (2006). A plant growth retardant, uniconazole, is a potent inhibitor of ABA catabolism in Arabidopsis. *Biosci. Biotechnol. Biochem.* 70, 1731–1739. doi: 10.1271/bbb.60077

Salas, R. A., Scott, R. C., Dayan, F. E., and Burgos, N. R. (2015). *EPSPS* gene amplification in glyphosate-resistant Italian ryegrass (*Lolium perenne* ssp. *multiflorum*) populations from Arkansas, USA. *J. Agric. Food Chem.* 63, 5885–5893. doi: 10.1021/acs.jafc.5b00018

Sasal, M. C., Demonte, L., Cislaghi, A., Gabioud, E. A., Oszust, J. D., Wilson, M. G., et al. (2015). Glyphosate loss by runoff and its relationship with phosphorus fertilization. *J. Agric. Food Chem.* 63, 4444–4448. doi: 10.1021/jf505533r

Serra, A. A., Couée, I., Renault, D., Gouesbet, G., and Sulmon, C. (2015). Metabolic profiling of *Lolium perenne* shows functional integration of metabolic responses under subtoxic conditions of chemical stress. *J. Exp. Bot.* 66, 1801–1816. doi: 10.1093/jxb/eru518

Serra, A. A., Nuttens, A., Larvor, V., Renault, D., Couée, I., Sulmon, C., et al. (2013). Low environmentally relevant levels of bioactive xenobiotics and associated degradation products cause cryptic perturbations of metabolism and molecular stress responses in *Arabidopsis thaliana*. *J. Exp. Bot.* 64, 2753–2766. doi: 10.1093/jxb/ert119

Sievers, F., Wilm, A., Dineen, D., Gibson, T. J., Karplus, K., Li, W., et al. (2011). Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. *Mol. Syst. Biol.* 7:539. doi: 10.1038/msb.2011.75

Singh, N., Megharaj, M., Kookana, R. S., Naidu, R., and Sethunathan, N. (2004). Atrazine and simazine degradation in Pennisetum rhizosphere. *Chemosphere* 56, 257–263. doi: 10.1016/j.chemosphere.2004.03.010

Smith, K., Simpson, R. J., Oram, R. N., Lowe, K. F., Kelly, K. B., Evans, P. M., et al. (1998). Seasonal variation in the herbage yield and nutritive value of perennial ryegrass cultivars with high or normal herbage water soluble carbohydrate concentrations grown in three contrasting Australian dairy environments. *Aust. J. Exp. Agric.* 38, 821–830. doi: 10.1071/EA98064

Soeno, K., Goda, H., Ishii, T., Ogura, T., Tachikawa, T., Sasaki, E., et al. (2010). Auxin biosynthesis inhibitors, identified by a genomics-based approach, provide insights into auxin biosynthesis. *Plant Cell Physiol.* 51, 524–536. doi: 10.1093/pcp/pcq032

Stamm, M. D., Enders, L. S., Donze-Reiner, T. J., Baxendale, F. P., Siegfried, B. D., and Heng-Moss, T. M. (2014). Transcriptional response of soybean to thiamethoxam seed treatment in the presence and absence of drought stress. *BMC Genomics* 15:1055. doi: 10.1186/1471-2164-15-1055

Steinrücken, H. C., and Amrhein, N. (1980). The herbicide glyphosate is a potent inhibitor of 5-enolpyruvylshikimic acid-3-phosphate synthase. *Biochem. Biophys. Res. Commun.* 94, 1207–1212. doi: 10.1016/0006-291X(80) 90547-1

Stiti, N., Missihoun, T. D., Kotchoni, S. O., Kirch, H. H., and Bartels, D. (2011). Aldehyde dehydrogenases in *Arabidopsis thaliana*: biochemical requirements, metabolic pathways, and functional analysis. *Front. Plant Sci.* 2:65. doi: 10.3389/fpls.2011.00065

Sulmon, C., Gouesbet, G., Couée, I., and El Amrani, A. (2004). Sugar-induced tolerance to atrazine in *Arabidopsis* seedlings: interacting effects of atrazine and soluble sugars on *psbA* mRNA and D1 protein levels. *Plant Sci.* 167, 913–923. doi: 10.1016/j.plantsci.2004.05.036

Sulmon, C., Gouesbet, G., El Amrani, A., and Couée, I. (2007). Involvement of the ethylene-signalling pathway in sugar-induced tolerance to the herbicide

atrazine in *Arabidopsis thaliana* seedlings. *J. Plant Physiol.* 164, 1083–1092. doi: 10.1016/j.jplph.2006.11.005

Sundin, L., Vanholme, R., Geerinck, J., Goeminne, G., Höfer, R., Kim, H., et al. (2014). Mutation of the inducible *Arabidopsis thaliana* Cytochrome P450 reductase 2 alters lignin composition and improves saccharification. *Plant Physiol.* 166, 1956–1971. doi: 10.1104/pp.114.245548

Teixeira, M. C., Duque, P., and Sá-Correia, I. (2007). Environmental genomics: mechanistic insights into toxicity of and resistance to the herbicide 2,4-D. *Trends Biotechnol.* 25, 363–370. doi: 10.1016/j.tibtech.2007.06.002

Valluru, R., and Van den Ende, W. (2011). Myo-inositol and beyond-emerging networks under stress. *Plant Sci.* 181, 387–400. doi: 10.1016/j.plantsci.2011.07.009

Velini, E. D., Alves, E., Godoy, M. C., Meschede, D. K., Souza, R. T., and Duke, S. O. (2008). Glyphosate applied at low doses can stimulate plant growth. *Pest Manag. Sci.* 64, 489–496. doi: 10.1002/ps.1562

Vigeland, M. D., Spannagl, M., Asp, T., Paina, C., Rudi, H., Rognli, O. A., et al. (2013). Evidence for adaptive evolution of low-temperature stress response genes in a Pooideae grass ancestor. *New Phytol.* 199, 1060–1068. doi: 10.1111/nph.12337

Voigt, C., Oster, U., Börnke, F., Jahns, P., Dietz, K. J., Leister, D., et al. (2010). In-depth analysis of the distinctive effects of norflurazon implies that tetrypyrrole biosynthesis, organellar gene expression and ABA cooperate in the GUN-type of plastid signalling. *Physiol. Plant.* 138, 503–519. doi: 10.1111/j.1399-3054.2009.01343.x

Walter, H., Consolaro, F., Gramatica, P., Scholze, M., and Altenburger, R. (2002). Mixture toxicity of priority pollutants at no observed effect concentrations (NOECs). *Ecotoxicology* 11, 299–310. doi: 10.1023/A:1020592802989

Wang, L., Hu, W., Sun, J., Liang, X., Yang, X., Wei, S., et al. (2015). Genome-wide analysis of *SnRK* gene family in *Brachypodium distachyon* and functional characterization of *BdSnRK2.9*. *Plant Sci.* 237, 33–45. doi: 10.1016/j.plantsci.2015.05.008

Ward, J. A., Ponnala, L., and Weber, C. A. (2012). Strategies for transcriptome analysis in nonmodel plants. *Am. J. Bot.* 99, 267–276. doi: 10.3732/ajb.1100334

Weisman, D., Alkio, M., and Colón-Carmona, A. (2010). Transcriptional responses to polycyclic aromatic hydrocarbon-induced stress in *Arabidopsis thaliana* reveal the involvement of hormone and defense signalling pathways. *BMC Plant Biol.* 10:59. doi: 10.1186/1471-2229-10-59

Yu, Q., An, L., and Li, W. (2014). The CBL-CIPK network mediates different signaling pathways in plants. *Plant Cell Rep.* 33, 203–214. doi: 10.1007/s00299-013-1507-1

Yu, Q., Han, H., Cawthray, G. R., Wang, S. F., and Powles, S. B. (2013). Enhanced rates of herbicide metabolism in low herbicide-dose selected resistant *Lolium rigidum*. *Plant Cell Environ.* 36, 818–827. doi: 10.1111/pce.12017

Zhang, Z. W., Feng, L. Y., Cheng, J., Tang, H., Xu, F., Zhu, F., et al. (2013). The roles of two transcription factors, ABI4 and CBFA, in ABA and plastid signalling and stress responses. *Plant Mol. Biol.* 83, 445–458. doi: 10.1007/s11103-013-0102-8

Zhang, Z. W., Yuan, S., Xu, F., Yang, H., Chen, Y. E., Yuan, M., et al. (2011). Mg-protoporphyrin, haem and sugar signals double cellular total RNAs against herbicide and high-light-derived oxidative stress. *Plant Cell Environ.* 34, 1031–1042. doi: 10.1111/j.1365-3040.2011.02302.x

Zhong, G. Y., and Burns, J. K. (2003). Profiling ethylene-regulated gene expression in *Arabidopsis thaliana* by microarray analysis. *Plant Mol. Biol.* 53, 117–131. doi: 10.1023/B:PLAN.0000009270.81977.ef

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Serra, Couée, Heijnen, Michon-Coudouel, Sulmon and Gouesbet. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.