Stacked genetically modified soybean harboring herbicide resistance and insecticide rCry1Ac shows strong defense and redox homeostasis disturbance after glyphosate-based herbicide application

World agricultural production of genetically modified (GM) products, in particular, the combination of different traits/genes in the same plant has been a trend over the last decade. There have been concerns raised over stacking multiple herbicide and insect-resistant transgenes that could result in fitness costs depending on the type and strength of selection pressures exerted by the environment. Here, we report the results of transcriptomic analysis comparing the effect of glyphosate-based herbicide (GBH) in the single-transgene versus stacked, herbicide-resistant soybean varieties on various biological processes, metabolic pathways, and key shikimic enzymes. Gene expression data showed that defense metabolism and redox homeostasis were equally modulated in single-transgene and stacked-variety samples. Carbon accumulation and energy metabolisms were distinct between the varieties and photosynthesis metabolism was found negatively affected in the single-transgene variety only. In the stacked variety, the shikimate pathway was modulated by the accumulation of transcripts from phenylalanine gene and other cascade genes. As expected, the expression of native EPSPS was upregulated in both varieties when herbicide was applied. On the other hand, transgenic EPSPS expression was down-regulated in both GM varieties upon herbicide application which cannot be explained. Glyphosate-based herbicides toxicity suggests its effects on carbon central metabolism and flux, redox metabolism, photosynthesis, and to hormone and defense response in plants. The observed unintended effects in GM herbicide-tolerant varieties unravel the deleterious effects previously observed on GM-tolerant varieties growth and production. The impact of GBH on shikimate and cascade pathways was observed in terms of both native and transgenic insensitive EPSPS modulation, alteration of jasmonic acid and lignin metabolism in both single-transgene and stacked variety. The energy metabolism and carbon flux were differently affected in these varieties. Oxidative stress, more specifically glutathione metabolism, induced by GBH, was also observed in this study. The stacked variety showed a more pronounced stress response (activation of specific stress defense proteins, Rboh, WRKY) and secondary compounds (β-glucosidase, isoflavone 7-O-methyltransferase). Omics profiling techniques, such as transcriptomics, can be considered tools to support risk assessment in detecting unintended effects due to the GBH application.

Background The combination of different traits or genes in genetically modified (GM) plants has rapidly emerged in worldwide crop production. In 2018, an increasing number of GM plants with stacked traits reached about 81 million hectares equivalent to 42% of the total 191.7 million hectares planted with transgenic crops worldwide [1]. The predominant trait, for both stacked and single-transgene crop varieties is herbicide resistance and it is estimated to remain so in the near future [2].
According to the current regulatory practice within the European Union (EU), stacked events are considered new GM organisms (GMO), requiring similar risk assessment procedures to those from single-transgene events [3], whereas in other countries, such as Brazil, stacked events are also considered new GMOs but require simplified risk assessments upon approval of single-transgene parental events [4].
Previous studies have shown that stacking herbicide and insect-resistant transgenes in Brassica sp. can result in fitness costs that are dependent on the type and strength of selection pressure, and could also impact plant communities through hitchhiking of unselected traits [5]. In that particular study, one of the tested selective pressure was the spray of GBH-, which has been shown to adversely affect plant uptake and transport of micronutrients (e.g. Mn, Fe, Cu, and Zn) and consequently, reduce disease resistance and plant growth [6,7].
Glyphosate, the active ingredient of Roundup herbicide works by inhibiting the enzyme 5-enolpyruvylshikimate-3-phosphate synthase (EPSPS), which catalyzes the penultimate step of the shikimate pathway leading to the conversion of shikimic acid to chorismate, the precursor for aromatic amino acids (tyrosine, phenylalanine, and tryptophan) and other secondary plant metabolites. Glyphosate competes with phosphoenolpyruvate (PEP), a substrate for the EPSPS enzyme, to form a very stable enzyme-herbicide complex that inhibits the product-formation reaction [8]. In glyphosate-tolerant crops, expression of transgenic CP4 EPSP synthase enables weed control by allowing post-emergent herbicide application [9], whereas in susceptible plants, EPSPS is inhibited by GBH thus causing a cascade of metabolic effects that are associated with glyphosate toxicity [10][11][12].
Despite the widespread use of GBH in global crop production, its precise mode(s)-of-action and potential adverse effects in GM plants remain unclear. Earlier studies with GM soybean varieties have shown a reduction in nodular activity, alteration in photosynthetic activity and a decrease in biomass when plants were sprayed with 4.8 kg a.i./ha of GBH [13]. The application of GBH in GM soybeans was also associated with phytotoxic responses, compromising essential biological processes and favoring the appearance of diseases [14]. Effects on secondary metabolism, oxidative, hormonal status, changes in photosynthesis were also observed for the same CP4 ESPSPS single cassette transgenic event [15][16][17][18].
The employment of transcriptomic and proteomic profiling techniques has supported the investigation of common metabolic responses to GBH applications in resistant plants. The accumulation of GBH in singletransgene varieties has enhanced cellular oxidation, possibly through mechanisms involving stimulation of the photorespiratory pathway [16]. In general, most of the GBH-induced genes are homologous to the known expression sequence tags-ESTs induced by abiotic stress factors [15].
Although previous studies provide insights into the genetic response of GM-resistant genotypes to GBH applications, an analysis of resistant genotypes conferring multiple transgenic traits is still lacking. This is partly because previous untargeted omics approaches were limited to the analysis of single-transgene, herbicide-resistant varieties [15,16,[19][20][21].
In order to gain an understanding about the impact of GBH applications in GM plants harboring two or more transgenic events, the current study analyzed several biological processes, metabolic pathways and key enzymes of the shikimate through transcriptomic profiling on the stacked INTACTA RR2 PRO soybean variety. We hypothesized that transgenic plants with a combination of transgenes respond differently to GBH exposure due to (1) the cost of expressing more than one heterologous protein and (2) synergistic and antagonistic interactions between transgenes and GBH target and cascade pathways. The data presented here provide new knowledge concerning the interactions between recombinant Cry1Ac (rCry1ac) and CP4-EPSPS transgene cassettes on the defense response and glutathione metabolism, the abundance of beta-glucosidase and oxidoreductase enzymes when GBH is applied. The experiment was conducted in a full-factorial random experiment in block design with two factors: soybean variety and herbicide treatment. Seedlings were grown in 14 L plastic pots filled with a substrate (1/3 clay soil; 1/3 cellulose residue and 1/3 poultry organic residue) with pH corrected to 6.0. The experiment was conducted in the greenhouse under two treatment conditions: herbicide spray application (treated group); and no herbicide application (control group). There were three plants per pot, and three pots per treatment, disposed in three random blocks and border protected. The experiment was watered daily. All plants were kept under the same conditions until they reached V2 stage, approximately 34 days after emergence, when the herbicide treatment was applied to a subset of plants (treated group). The treatment was conducted through spray application using GBH formulation Roundup Transorb ® (Monsanto do Brasil S.A.) (Potassium salt of N-(phosphonomethyl) glycine (588 g/L; N-(phosphonomethyl) glycine acid (480 g/L; Other ingredients (820 g/L). The maximum dosage informed by the leaflet (4.5 L/ha; 2.2 kg a. i./ha) for soybean crops was applied [22]. To minimize spray contamination and drift, GBH was applied to treated samples outside the greenhouse. After 8 h, the pots were returned to the greenhouse to collect leaf material [12]. For each treatment, three leaves of the fourth trifoliate were collected from the three plants present in each pot, immediately frozen in liquid nitrogen and stored at − 80 °C. Hence, each treatment contained three biological pools of three different plants, which were used for the transcriptomic analysis (Fig. 1).

Library construction, sequencing and mapping
Total RNA was isolated using RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The purity of RNA was measured using Nanodrop 2000 (Thermo Scientific), and the integrity of RNA was determined by agarose gel electrophoresis. Samples of 1-2 μg mRNA were used for library construction and sequencing. Libraries were constructed by Novogene Corporation (Sacramento, CA). Briefly, mRNA was enriched from total RNA using oligo (dT) beads. The mRNA was then randomly fragmented and cDNA was synthesized using random hexamers. Library construction consisted of terminal repair, A-tailing, ligation of sequencing adapters and size selection and PCR enrichment. Library concentration was determined and normalized to 1 ng/μl. Libraries were sequenced on Illumina HiSeq platform (PE150). Raw reads are filtered to avoid reads containing adapters or low quality. The Hisat2 v2.1.0-beta algorithm was used to map-filtered sequenced reads into Glycine max reference genome (Glycine_max_v2.1) [23]. The quantity of total mapped reads and its percentage of clean reads were calculated and they are presented in Table 1. The transcriptome raw sequencing data from this study have been submitted through the Sequence Read Archive (SRA) on the NCBI database (http://www.ncbi.nlm.nih.gov/) under the BioProject number PRJNA625648.
Mapped transcripts were annotated and counted generating read count values in terms of fragments per kilobase of transcript per million mapped reads (FPKM) by using HTSeq v0.6.1 software (union mode). Differentially expressed genes (DEGs) between the different treatments were determined by using read counts from gene expression level analysis with DESeq v1.10.1 software following three main steps: read counts normalization; p value estimation (negative binomial distribution); and false discovery rate (FDR) value estimation [24]. Only genes with FDR adjusted p value (q value) < 0.05 after Benjamini-Hochberg correction for multiple-testing [25] were considered as significant DEGs.

Functional annotation of DEGs and bioinformatics
DEGs were annotated for gene ontology (GO) terms using Blast2GO v2.5 and categorized into Molecular Function (MF), Cellular Component (CC), and Biological Process (BP) categories [26]. DEGs were also annotated using KEGG pathway enrichment analysis (Kyoto Encyclopedia of Genes and Genomes) aiming at identifying significantly enriched metabolic pathways or signal transduction pathways affected by the GBH treatment. Pathways with q-value < 0.05 were significantly enriched. A heatmap clustering analysis of the log10(FPKM + 1) values was conducted using pheatmap library in R environment aiming to find gene expression patterns across the different treatments [27,28].

Epsps transcript quantification by qRT-PCR
qRT-PCR analysis was performed aiming at quantifying the expression of native and transgenic epsps transcripts in both single-transgene and stacked soybean varieties. DNA primers were designed based on differences in coding sequence of each of the two targets: epsps transgene in MON89788-1 and MON4032-6 and native epsps (PrimerQuest; Integrated DNA Technologies Inc., Skokie, IL, USA) (Additional file 1). The cDNA was synthesized using 100 ng of total RNA from each biological replicate (Superscript VILO cDNA synthesis kit; Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Real-time PCR reactions were performed using tenfold diluted cDNA product and set up using the Power SYBR green PCR master mix (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's directions and run on a StepOnePlus Real-Time PCR system (Applied Biosystems). Thermocycling conditions were: 50 °C for 2 min; 95 °C for 10 min; and 40 cycles at 95 °C for 15 s, 60 °C for 1 min, with melt curve set at 95 °C for 15 s, 60 °C for 1 min, 95 °C for 30 s, and 60 °C for 1 min. Five housekeeping genes were tested for normalization calculation: B-actin [29], F-box [30], 60S [31], 18SrRNA [32] and Skip16 [33]. Normalization of cycle threshold (Ct) values was performed using the 60S ribosomal protein housekeeping gene as it showed low variability (< 15%) and high-efficiency results (98,59%; R 2 = 0.99). Relative expression levels (Rq) were obtained using ΔΔCt method (StepOne Plus Software v2.3; Applied Biosystems). Statistical significance was assessed through t test (p < 0.05) [34]. This experiment was performed twice (new seedlings grown) in order to confirm the results. In the first experiment, samples from the RNA-Seq experiment were used (three biological replicates consisting of a pool of three seedlings in each of the three randomized blocks in a total of 27 seedlings analyzed per treatment). A second experiment consisted of new seedlings (six biological replicates consisting of a pool of three seedlings in each of the three randomized blocks in a total of 54 seedlings analyzed per treatment).

Validation of gene expression patterns by RT-qPCR
In order to further verify the expression profiles generated by Illumina sequencing analyses, 16 transcripts were selected for RT-qPCR. The same PCR conditions were applied from the previous epsps quantification experiment, including the housekeeping gene. Transcripts were selected based on their contribution to biochemical pathways most affected by GBH treatment (i.e., Photosynthesis; Plant hormone signaling; Glutathione metabolism; Phenylpropanoid biosynthesis; Flavonoid and Isoflavonoid biosynthesis; Monoterpenoid biosynthesis; Protein processing in endoplasmic reticulum; Phagosome; and Propanoate metabolism). The relative quantity (Rq) obtained by RT-qPCR for selected transcripts is shown in Additional file 2. Rq values were calculated as the relative expression in treated plants compared to control plants (Rq = 1). Log2FC values obtained from RNA-seq analysis for the same transcripts are also presented in Additional file 2 for comparative purposes. All transcripts analyzed by RT-qPCR showed highly similar expression trends compared to RNA-Seq data. These results validate the sequencing data produced in this study, as well as indicate that such data are accurate and reliable. A total of 1425 (1024 up-regulated; 401 down-regulated) and 547 transcripts (522 up-regulated; 25 downregulated) were identified as DEGs in response to herbicide treatment for the single-transgene and stacked varieties, respectively (Fig. 2). A complete list of all DEGs and corresponding p-adjusted values are available in Additional file 3.

Transcriptome assembly and gene expression
We have further explored the data by running a hierarchical clustering analysis of DEGs aiming to find genes with similar expression patterns across the different varieties and treatments. The heatmap showed that gene expression data clustered according to treatment, meaning the factor 'herbicide' resulted in a major effect compared to the effect of the genetic background-'variety' factor ( Fig. 3). Although no precise information regarding the genetic background of the two commercial hybrids used were obtained, such result is expected once the soybean genetic diversity in Brazil is relatively low. For example, most of the Brazilian soybean germplasm is derived from four main genotypes (CNS, S-100, Roanoke and Tokyo), which contributed to more than a half of the . Data for genes that were not classified as differentially expressed are plotted in blue. In green and red, we plotted data for genes that are differentially expressed after application with glyphosate-based herbicide (Bonferroni corrected p-value ≤ 0.05) with an absolute log2 fold change (|FC|) greater than 1.5 genetic base of all commercial cultivars released in Brazil [35,36].

Gene Ontology annotation of DEGs
In the single-transgene variety, herbicide application resulted in 16 significant BP terms being up-regulated, with protein phosphorylation at highest modulation. A total of six BP terms were annotated as significantly down-regulated and the two most enriched terms were photosynthesis and response to auxin (Fig. 4a). Similarly, the stacked variety showed 10 BP terms up-regulated, all related to protein metabolism. No significant BP term was found to be down-regulated (Fig. 4b). Interestingly, distinct DEGs for MF terms were observed for singletransgene and stacked varieties in response to herbicide application. The single-transgene variety showed protein kinase activity term as the most up-regulated term (1/13 terms), and copper ion binding as the only significant down-regulated MF term. The stacked variety showed a total of eight MF terms that were significantly up-regulated, in which the most enriched terms were related to catalytic activity (i.e. oxidoreductase and fatty acid) and nucleic acid binding. No terms for down-regulated DEGs were significantly annotated for the stacked variety. As for the Cellular Component domain, the single-transgene variety showed endoplasmic reticulum term annotated as up-regulated, while eight significant down-regulated terms were related to the Photosystem II. The stacked variety showed nuclear chromosome and extracellular matrix CC terms being up-regulated, and no terms were annotated as significantly down-regulated.

Pathways with differentially expressed genes
Pathway enrichment analysis of DEGs was performed in order to test whether GBH induces metabolic changes in stacked-transgene soybean varieties compared to singletransgene varieties at official dosage recommendations (2.2 kg a.i. ha). The list of all KEGG pathways and corresponding annotated genes are available in Additional file 4.
The most affected pathways in the single-transgene variety were Protein processing in endoplasmic If there are less than 20 pathways, all the pathways will be presented. X-axis represents the name of the pathway and the Y-axis represents the Rich Factor. Rich factor is the ratio of the DEG number to the background number in a certain pathway. The size stands for the number of difference genes and the color stands for different q-values reticulum and Protein export, Plant-pathogen interaction, Biosynthesis of secondary metabolites, Monoterpenoid biosynthesis, Carbon metabolism and Biosynthesis of amino acids (Fig. 5a). In the stacked variety, Biosynthesis of secondary metabolites, Plant-pathogeninteraction, Cyanoamino acid metabolism, Protein processing in endoplasmic reticulum, Phenylalanine metabolism and Cysteine and methionine metabolism (Fig. 5b).

Impact of GBH in the shikimate pathway and direct cascade effects
Despite the widespread use of GBH in agriculture, major questions remain on how its exposure affects cell metabolism and physiology in glyphosate-resistant plants and if there are antagonistic or synergistic effects in stackedtransgene varieties. In order to address these questions, we have profiled transcriptomic changes after GBH treatment and characterized the interactions between the shikimate pathway and other unsupervised metabolic pathways. We were also interested in the direct effects of GBH in native and transgenic EPSPS expression.
Our results showed that native EPSPS expression showed a 1.5-fold change up-regulation in both single-transgene and stacked events, which suggests that reduced levels of AAAs in response to EPSPS inhibition may act as a signal to induce the expression of the shikimate pathway genes and restore the carbon flux through the pathway in plants [37]. On the other hand, insensitive transgenic cp4-epsps showed a decrease in transcript accumulation also in both transgenic varieties. Although stably integrated into the genome, variable and non-directional levels of CP4 EPSPS have been observed linked with other factors, such as the genetic background, trait stacking, growing region or season [38]. But the extent to which detection protocols could differentiate both versions of EPSPS is unclear. We were careful in repeating the full epsps experiment with new seedlings and new herbicide application in order to confirm these results. The underlying mechanism for the decreased While EPSPS enzymes have been equally modulated in both single-transgene and stacked varieties, differences found in downstream metabolism between the varieties suggest synergistic and antagonist effects of GBH when transgenes are stacked. In the stacked variety, the modulation of other amino acids metabolism (cysteine and methionine) and secondary metabolites metabolism, such as flavonoids and glutathione metabolism, as well as the jasmonic acid metabolism, were most prominent (Fig. 7). Jasmonic acid derives from fatty acid biosynthesis and increased levels have been also observed after GBH and drought stress application in NK603 herbicideresistant GM maize [20]. GBH was also shown to affect other hormones, such as ethylene [39] and abscisic acid [12].
Intermediates in the shikimate pathway are used for the synthesis of proteins and that in plants also serve as precursors of numerous natural products, such as pigments, alkaloids, hormones, and cell wall components. [40]. Such plant natural compounds play crucial roles in plant growth, development, reproduction, defense, and environmental responses [41]. For example, phenylalanine is a common precursor of numerous phenolic compounds, including lignin, which corresponds to the highest carbon flux [37]. We observed a 3.5-fold change increase of PAL, an enzyme involved in lignin biosynthesis, in the stacked variety. Differently, Zobiole et al. [42] found negative correlation between levels of lignin in single-transgene GM soy and increasing rates of GBH applications. Previous studies also reported higher susceptibility levels to diseases after application of GBH in transgenic varieties which were associated with changes in lignin contend and, consequently, with morphological and functional quality of the plant defense organs [6,43,44].
Toxicity of GBH applications is a known and desired outcome in susceptible plants, whereas in GM-tolerant plants, the modulation of shikimate genes is not expected as they contain insensitive and constitutively expressed CP4-EPSPS. Our results show that GM single-transgene and stacked varieties have strong modulation of the shikimate pathway, including EPSPS versions, which also affected numerous compounds derived from shikimate precursor chorismate molecule.

Changes in central carbon metabolism and carbon flux
The constitutive expression of transgenes controlled by strong viral promoters, such as P35S from the cauliflower mosaic virus, has been always a concern due to potential energy cost and carbon flux in plant cells. In this paper, we applied GBH, an inhibitor of the enzyme EPSPS, at recommended concentrations in order to investigate its Fig. 7 Schematic representation of modulated metabolic pathways in stacked-transgene soybean variety treated with GBH. Legend: The shikimate pathways support the formation of numerous natural products in plants. It produces chorismate, a common precursor for the tryptophan (Trp) pathway, the phenylalanine/tyrosine (Phe/Tyr) pathways, and the pathways leading to folate, phylloquinone, and salicylate. Trp, Phe, and Tyr are further converted to a diverse array of plant natural products that play crucial roles in plant physiology. The gray boxes show metabolic pathways altered in this study. The blank boxes are metabolic pathways found altered in the literature. Filled arrows demonstrate indirect relationships between the metabolic routes and the shikimate pathway, whereas dashed lines show a direct link impact in downstream metabolism in GM-tolerant plants containing two or more inserts. Under normal growth conditions, more than 30% of plant-fixed carbon flows through the shikimate pathway [37,45,46], whereas under stress, plants mobilize their carbon stocks to transform energy and resist to harmful effects.
In susceptible plants, inhibition of EPSPS reduces the levels of aromatic amino acids and their downstream products which act as a signal to induce the expression of the shikimate pathway genes and restore the carbon flux through the pathway. However, the effects of GBH in the central carbon metabolism of GM-tolerant plants have also been previously observed.
In our study, we observed a clear difference in carbon metabolism between the single-transgene and the stacked variety. While single-transgene variety showed modulation in carbon fixation and glycolysis metabolism, the stacked variety exhibited changes in starch and sucrose metabolism. Our data suggest that single-transgene varieties enhance the cytosolic glycolytic network to provide metabolic flexibility that facilitates plant acclimation to herbicide stress. Modulation of phosphoenolpyruvate carboxykinase that promotes reversible protein phosphorylation of major importance in controlling legume nodule carbon metabolism and related metabolite transport was observed. In addition, other enzymes involved in parallel reactions were also altered; fructose bisphosphate aldolase and malate dehydrogenase genes were also found up-regulated. These changes are in agreement with previous studies that have observed changes in GM soybean nodulation after GBH applications [13,47].
On the other hand, increased levels of trehalosephosphate transcripts (3.6 log2FC) were observed in the stacked variety. Sucrose and starch balance is directly related to optimization of growth rates [48,49]. Trehalose (α-d-glucopyranosyl-1,1-α-d-glucopyranoside) works as an osmolyte, storage reserve, transport sugar, and stress protectant [50,51]; and it is also involved in growth and development metabolism [52] with clear links to abscisic acid and auxin signaling [53]. Increased levels of trehalose have been observed in response to osmotic stress [54] as well as to dehydration stress tolerance [54,55]. Most plants accumulate substantial starch reserves in their leaves to provide carbon and energy for maintenance and growth [56,57]. Therefore, the accumulation of soluble sugars, such as trehalose, is suggested to be a protective mechanism under oxidative stress conditions [58,59].

Altered cellular redox homeostasis
Exposure to GBH is directly linked to accumulation of antioxidant enzymes, indicating an oxidative stress [60]. Glutathione (GSH) is a key molecule in the antioxidant network in plants, acting to control reactive oxygen species (ROS) accumulation and facilitating cellular redox homeostasis especially under stress conditions [61]. In fact, GSH plays an important role in herbicide detoxification via the glutathione S-transferase (GST) system [62]. We found evidence for cellular detoxification response through significant up-regulation of GST in both transgenic varieties (single-transgene variety: average log2FC = 3.1; Stacked: average log2FC = 3.5).
Contrarily, other genes encoding important enzymes related to glutathione metabolism showed to be differently affected in the single-transgene and stacked varieties, revealing that both genotypes may respond differently to oxidative stress. For instance, glucose 6-phosphate dehydrogenase (G6PDH), an enzyme participating in the first two reactions of oxidative pentose phosphate pathway, was significantly down-regulated in the singletransgene variety. Reduced levels of G6PDH is related to glutathione depletion and subsequent high oxidative stress in the cell [63]. Consequently, reduced glutathione (GSH) is required to combat oxidative stress and maintain the normal reduced state in the cell, a phenomenon known as the redox homeostasis [16,60,61]. Oxidized glutathione (GSSG) is reduced to GSH by NADPH generated by G6PDH in the pentose phosphate pathway [64]. Complete depletion of glutathione in its reduced form (GSH), or the production of GSSG from GSH, with concomitant accumulation of formaldehyde have already been reported as signs of undergoing oxidative stress in single-transgene soybean event as compared to its non-GM isogenic line [19,65].
In the stacked variety, although G6PDH gene expression has not been significantly affected, herbicide treatment up-regulated the expression of 6-phosphogluconate dehydrogenase (6PGDH) gene (log2FC = 1.25). 6PGDH, a second enzyme participating in the OPPP, catalyzes the NADP-dependent oxidative decarboxylation of 6-phosphogluconate generating NADPH and ribulose-5-phosphate, a precursor for the synthesis of nucleotides and nucleic acids [66]. We hypothesize that the production of such reducing equivalents is being used in further reductive reactions in stacked plants, such as keeping GSH in its reduced form, aiming at maintaining the cell redox homeostasis.
Our results also showed protein processing in endoplasmic reticulum (ER) as one of the most up-regulated pathways in both, single-transgene and stacked varieties when GBH is applied. Glutathione homeostasis in response to oxidative stress has been also described as active in the ER [67]. A diverse range of genes encoding important molecular chaperones guiding secretory folding proteins, as well as ubiquitin-proteasomes responsible for exporting and degradation of misfolded proteins, were shown to be significantly up-regulated in the presence of GBH. Since glutathione is oxidized, transport proteins must export GSSG from the ER to the cytosol aiming to reach an ideal glutathione homeostasis [67]. Conversely, the stacked variety showed evidence of oxidative stress responses due to the up-regulation of cytosolic glutathione genes (GST log2FC = 3.5; 6PGDH log2FC = 1.25), while only genes encoding ERAD enzymes were significantly up-regulated in ER. Vivancos et al. [16] have also found effects of herbicide on cellular redox homeostasis of single-transgene, GBHresistant soybean varieties. More specifically, the authors reported that the accumulation of high levels of glyphosate in GM plants-enhanced cellular oxidation, possibly through mechanisms involving increasing of photorespiratory pathway [16]. Moreover, a recent integrative in silico model of C1 metabolism in single-transgene, GBHresistant GM soybean predicted complete depletion of glutathione and accumulation of formaldehyde as a result of oxidative stress compared to its non-GM counterpart [19]. According to our findings, single-transgene and stacked GM soybean showed oxidative stress at different levels and cellular components.

Photosynthesis imbalance
Photosynthesis efficiency and inhibition of chlorophyll function has been observed as a side-effect from GBH applications in both susceptible and GM-tolerant plants [16,18,68]. In other words, GBH seems to impact photosynthesis as a side-effect of glyphosate and its byproducts and/or adjuvants whether or not insensitive epsps sequences are present in the genome. In our study, the single-transgene variety showed a decrease in the light-harvesting chlorophyll A and B content (complex I of class LhcA 2,3 and 4 with four genes involved, and the complex II of class LhcB 1,2,3 and 6 with nine genes involved). These findings are supported by Li et al. [69], which also observed a decline in the content of chlorophyll A and B in GM and conventional soybean varieties under GBH treatment [69]. In addition, we found two genes down-regulated related to putative ferredoxin enzymes. The amount of ferridoxin is also decreased in tobacco under various stresses, including those from herbicide treatment [70]. Iquebal et al. [71] observed that genes involved in the photosynthetic pathway were deregulated after exposure to herbicides in resistant chickpea variety [71]. In Lolium perenne sensitive plants, chlorophyll fluorescence was also affected by glyphosate [72].

Defense and environmental responses
Defense imposes a substantial demand for resources that can negatively impact growth and diminish the overall set of energy reserves and/or promote resource diversion for growth, defense, and reduction of photosynthesis [73]. Previous transcriptome studies using microarray technique to investigate the metabolic impact of GBH treatment in susceptible and resistant soybean, arabidopsis and brassica showed that most affected pathways are involved in defense metabolism [12,15,74]. In this study, the single-transgene and stacked varieties showed up-regulation of calcium-related pathways, which are essential to coordinate a rapid adaptive response in several species [75,76]. Calmodulin protein families were also altered in both single-transgene variety (two altered genes, average 2.10 log2FC) and stacked (five altered genes, average 3.3 log2FC) varieties. Previous studies with GBH application in sensitive soybean also observed changes in calcium-related genes regardless of herbicide concentrations and the collection time after application (4 and 24 h) [12,29]. The Ca2+/CaM complex play key roles in plant metabolism as it is the main signal transduction pathway involved in turgor regulation and with an impact in drought tolerance [77].
Rapid recognition of injuries by cellular signal transduction pathways occurs through various signaling molecules, including calcium, protein phosphorylation and ROS, which are well-known triggers of stress resistance in plants [78]. Herbicides are considered abiotic stressors that can disrupt the balance between the production and elimination of ROS [79]. There is a close relationship between calcium-dependent ROS production and a specific group of genes. For example, the respiratory burst oxidase homolog (Rboh) gene family. Activation of this group occurs after the recognition of pathogens and a variety of other processes [80,81]. We observed strong up-regulation of the Rboh group (3.58 log2FC) in the stacked variety. Such oxidases have been reported as key factors in activating innate and mobilized immunity during oxidative stress damage [82].
Another example of defense regulatory circuit was the identification of WRKY transcription factors, which are connected to phosphorylation events of mitogen-activated protein kinase (MAPK) in response to pathogen recognition with the accumulation of Rboh protein [81] Strict regulation and fine-tuning of WRKY proteins are directly linked to plant stress signaling responses [83,84], such as saline stress [85], drought [86,87] and heat stress [88]. We observed up-regulation of WRKY genes in both varieties, with higher expression and number of genes in the stacked variety (five genes with an average of 2.5 fold change).
In addition, pathogenesis-related proteins (PR), known as an indispensable component of innate immune responses in plants under biotic or abiotic stress conditions, were also altered in this study. In the single-transgene variety, we find one gene PR1, up-regulated with a 2.9 log2FC. These proteins are also involved in hypersensitive response or systemic acquired resistance against a variety of plant infections [89] and an important response mechanism to multiple stresses [90]. PR proteins are considered the signature genes of salicylic acid and jasmonic acid pathways in many crop plants [90][91][92][93]. Furthermore, Hsp90 gene families were found up-regulated in the both single-transgene and stacked varieties. In soybean varieties, Hsp90 gene was found induced by heat, salt, and osmotic stresses [94]. Strikingly, the stacked variety up-regulated seven β-glucosidase-related genes with an average of 5 log2FC. We also found a regulated isoflavone 7-O-methyltransferase gene with 8.5 log2FC. β-glucosidases enable the enzymatic removal of a protecting glucose group, thus providing plants with an immediate chemical defense against protruding herbivores and pathogens [95].

Relevance to risk assessment of stacked GM crops
Worldwide, a growing number of GM crops with stacked transgenic traits are being developed to confer resistance to herbicide active ingredients and some insect species. For most varieties, the single-transgene events might never reach market and pre-market risk assessment. Therefore, an assessment of the risks of the actual GMO to be released in the environment should consider combinatorial and cumulative effects derived from stacking transgenes into single organism.
Omics profiling analysis can contribute to the identification of combinatorial effects that may occur due to interactions among the proteins and metabolites produced by the transgenes or endogenous genes of a stacked GM plant. In addition, interactions between the stacked transgenes or their products, or interactions among the physiological pathways in which the transgenes are involved, taking into account the possibility that these interactions could result in potentially harmful substances, such as anti-nutritional factors, some of which may persist or accumulate in the environment.
Stacked GM plants can be produced through different approaches. In addition to the cross-breeding of two GM plants, multiple traits can be also achieved by the natural cross of transgenic lines that have been found in crop field boundaries [96,97], such as feral transgenic canola outside of cultivation [98,99].
Accordingly, it is reasonable to anticipate future occurrence of stacked traits within ruderal and wild populations. Despite the potential for the formation of feral populations with multiple transgenes, we have little understanding of how these traits could migrate, evolve or influence native and naturalized plant communities. Thus, such profiling studies could generate useful information to assist risk assessment of stacked GM crops and potential feral populations.
Our study on stacked-transgene soybean variety showed GBH effects on shikimate genes, carbon metabolism and flux, photosynthesis, oxidative events and defense response. Whereas GBH effects in singletransgene plants have been reported [100][101][102] our data suggest that GBH affects stacked-transgene plants at a higher extend than its single-transgene near-isogenic comparator. In addition, GBH adverse effects in GM tolerant plants are widely variable depending on species and cultivar and herbicide regime. Therefore, it is intrinsically important to elucidate the GBH effects on physiological processes related to metabolic disturbances in order to better understand the glyphosate-herbicidal mechanism and its possible unintended effects on commercialized transgenic varieties.

Conclusions
It is well known that glyphosate kills undesired plants by inhibiting EPSPS synthase enzyme, thus blocking the synthesis of aromatic amino acids. However, glyphosatebased herbicides have shown to promote several indirect effects on plant physiology which may also explain its herbicidal effects. Glyphosate-based herbicides toxicity suggests its effects on carbon central metabolism and flux, redox metabolism, photosynthesis, and to the plant's hormone and defense response. Most relevant, such unintended effects are also present in GM herbicide-tolerant varieties even when they do not lead to plant's death. The alteration of these cellular processes unravels the deleterious effects previously observed on GM tolerant varieties growth and production. The impact of GBH on shikimate and cascade pathways was observed in terms of both native and transgenic insensitive EPSPS modulation, alteration of jasmonic acid and lignin metabolism in both single-transgene and stacked variety. Whereas the energy metabolism and carbon flux were differently affected in these varieties. In the stacked variety, trehalose levels were altered up to sevenfold increase. Oxidative stress, more specifically glutathione metabolism, induced by GBH, was observed in this study. Redox imbalance is known to severely damage the cell integrity and can negatively interfere with photosynthetic processes, for example by decreasing the chlorophyll content, photochemical efficiency, and C metabolism, leading to reduction in plant growth. We found Ca2 + signaling responses and several up-regulated molecular chaperones in both varieties. However, distinct stress responses were also observed. The stacked variety showed a more pronounced stress response (activation of specific stress defense proteins (Rboh, WRKY) and secondary compounds (β-glucosidase, isoflavone 7-O-methyltransferase). Unintended effects of GBH applications in GM tolerant varieties as well as the differences in the variety's response show the relevance in elucidating the GBH effects on physiological processes as means to establish its safety.