Transcriptome analysis reveals different responses and tolerance mechanisms of EPSPS and GAT genes in transgenic soybeans

Background: Glyphosate is a broad-spectrum and non-selective systemic herbicide. Introduction of glyphosate tolerance genes like EPSPS or detoxication genes like GAT could confer glyphosate tolerance to plants. Our previous study revealed that co-expression of EPSPS and GAT genes conferred higher glyphosate tolerance with no “yellow ashing” phenomenon. However, the tolerance mechanisms response to glyphosate at transcriptional level was still unknown. Methods and Results: To investigate glyphosate tolerance mechanisms in different GM soybeans, RNA-seq analysis was conducted using four soybean genotypes, including two non-transgenic (NT) soybeans ZH10, MD12, and two GM soybeans HJ698 and ZH10-6. A total of 90.72 Gb clean reads were generated and differentially expressed genes (DEGs) were identied in these soybeans before and after glyphosate treatments. The similar response of glyphosate in NT soybeans and the different effect of glyphosate to two GM soybeans were identied. The number of DEGs involved in shikimate biosynthetic pathway and herbicide targeted cross-pathways was continuously increased in NT soybeans, and slightly increased in HJ698 as the increasing of treat time. However, it was only signicantly increased at 12hpt but tended to recover to levels of 0 hpt at 72 hpt in ZH10-6, which can explain the higher glyphosate tolerance of ZH10-6. Conclusions: All these results suggested that GAT and EPSPS genes associated play a crucial role in response to glyphosate, and the former might work at the early stage of glyphosate exposure while the latter might be activated after the uptake of glyphosate uptake. These ndings will provide valuable sights for understanding of the molecular basis underlying glyphosate tolerance. were selected for transcriptome analysis. Transcriptional changes of these soybean genotypes contrasting in their tolerance to glyphosate at three time points were revealed by RNA-seq and thousands of differentially expressed genes were identied. The objective of this study was to investigate the gene expression proles of these soybeans, and to explore the possible mechanism in response to glyphosate. The results will provide not only a basis for understanding the mechanism in response to glyphosate, but also valuable clues for further glyphosate research in crops. experiments.


Introduction
Soybean (Glycine max L.) is the most frequently cultivated grain legume crops all over the world [1].With the development of transgenic technology, commercially available genetically modi ed (GM) soybeans are mostly tolerant to glyphosate or glufosinate.Glyphosate is the most widely used broad-spectrum, non-selective herbicides.The plant endogenous 5-enolypyruvylshikimate-3-phosphate synthase (EPSPS) activity is highly inhibited by glyphosate, leading to the accumulation of shikimate and the de ciency of chorismate in the shikimate pathway [2].However, the chorismate is the indispensable precursor for biosynthesis of aromatic amino acids in plants, including phenylalanine, tyrosine and tryptophan [2,3].
To date, majority of glyphosate-tolerant GM soybeans are derived from soybean event GTS 40-3-2, which contains one functional insert expressing CP4-EPSPS gene (http://www.isaaa.org).CP4-EPSPS encodes 5-enolypyruvylshikimate-3-phosphate synthase (EPSPS) whose low a nity for glyphosate and high catalytic activity contributes to glyphosate tolerance in plants [4].Some studies reported that overexpression of exogenous EPSPS conferred glyphosate tolerance in plants, but the applied glyphosate was found to accumulate in non-transgenic (NT) and GM plants [5][6][7][8].In turn, its accumulation led to the decrease of crop yield by interfering with photosynthesis, chlorophyll content, shikimate content and nutrient accumulation [9].The glyphosate tolerance does not mean the reducing of glyphosate translocation or accumulation in GM plants, and the main concern is the toxic injury caused by glyphosate accumulation.Previous studies also showed that the glyphosate accumulation led to the decrease of chlorophyll content, accumulation of primary phytotoxic metabolite or the formation of insoluble glyphosate-metal complexes [9][10][11].Interestingly, the typical symptom of glyphosate application, which is known as "yellow ashing", was observed in not only NT soybean but also GM soybean only expressed with EPSPS gene [9].With the identi cation of more glyphosate-related genes, the glyphosate N-acetyltransferase gene (GAT) or glyphosate oxidoreductase gene (GOX) was used in detoxi cation of glyphosate [12][13][14].Co-expression of GAT and EPSPS genes conferred higher tolerance to glyphosate in tobacco, cotton and soybean, and the "yellow ashing" symptom was not found after glyphosate application [8, 15,16].Moreover, the signi cant changes of shikimate acid content and chlorophyll content were not detected in GM soybean expressed with EPSPS and GAT genes, suggesting that the shikimate pathway could not be inhibited [16].
Over the past years, some reports have focused on the glyphosate-response mechanism and gene expression pro ling pattern with glyphosate treatment in NT soybean and GM soybeans expressed with CP4-EPSPS [17,18].However, little is known about how plants response to glyphosate when glyphosate exposure, particularly the transcriptional changes in plants expressed glyphosate tolerance EPSPS and detoxi cation gene GAT.In this study, two NT soybeans (ZH10 and MD12) and two GM glyphosatetolerant soybeans (ZH10-6 and HJ698) were selected for transcriptome analysis.Transcriptional changes of these soybean genotypes contrasting in their tolerance to glyphosate at three time points were revealed by RNA-seq and thousands of differentially expressed genes were identi ed.The objective of this study was to investigate the gene expression pro les of these soybeans, and to explore the possible mechanism in response to glyphosate.The results will provide not only a basis for understanding the mechanism in response to glyphosate, but also valuable clues for further glyphosate research in crops.

Plant materials and growth conditions
Glyphosate-tolerant soybeans ZH10-6, HJ698 and NT soybean cultivars Zhonghuang10 (ZH10), Mengdou12 (MD12) were used for glyphosate treatment and transcriptome sequencing.Among them, ZH10-6 expressed with glyphosate tolerance gene G2-EPSPS and glyphosate-degrading gene GAT was developed by Agrobacterium-mediated transformation of ZH10 [16].HJ698 is a backcross line derived from MD12 and GTS40-3-2.Soybean seeds were planted in 10 cm plastic pots containing mixture of vermiculite, peat, nutrient soil with the ratio of 1:2:4 and grown in the controlled chamber under a 16h/8h photoperiod with 120 μmol/m 2 •s light intensity.Day and light temperatures were 26℃ and 24℃, respectively.Three plants per pot were reserved after seedling emergence.

Glyphosate treatments
Commercially formulated isopropylamine salt of glyphosate with the rate of 300 g a.e./L (RoundupÒ, Monsanto Company) was used for glyphosate treatment.The labeled rate for glyphosate application varies from 600 to 1200 g a.e./ha [39].When the rst trifoliolate leaves of seedlings were fully expanded, soybean plants were spraying with RoundupÒ at the rate of 0.18 g a.e./m 2 as reported previously [16], which is a relatively high dosage represent the "worst-case scenario" to promote injury.The phenotype of whole seedlings were observed before (0 hpt) and 12, 24, 72, and 120 hours post treatment (12,24,72 and 120 hpt).At the same time, leaves from ten plants of these four soybean genotypes were collected at 0hpt, 12 hph and 72 hpt.All samples were immediately frozen in liquid nitrogen and stored at -80℃ until RNA extraction.

RNA extraction and sequencing
Total RNA was extracted from leaves with RNA extraction kit (TRIzol reagent, Invitrogen Inc) following the manufacturer's instructions.The quality and quantity of RNA were assessed with agarose gels electrophoresis and a Nanodrop spectrophotometer (Thermo Scienti c, USA).A total of 14 libraries were constructed using reagents available in the RNA-Seq Sample Preparation Kit (Illumina, USA).All these libraries were sequenced and paired-end reads were generated using the Illumina HiSeq 2500 platform.

RNA-Seq data ltering and analysis
Quality control of raw reads was performed using BMKCloud software (https://www.biocloud.net/)and high quality clean reads were obtained through removal of adaptor sequences, ambiguous reads with 'N' bases and low quality reads in which >10% of bases had a Q-value less than 20.Then all clean reads were aligned to the soybean reference genome of cv.Williams 82 (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Gmax) using TopHat software (http://ccb.jhu.edu/software/tophat/index.shtml).After passing the process of the alignment, the duplicate reads were removed and the remaining reads were used for further analysis.

Differential gene expression analysis
Gene expression analysis was performed by Cu ink 2.2.1 (http://cole-trapnelllab.github.io/cuinks/releases/v2.2.1/).Fragments per Kilobase of transcript per Million fragments mapped (FPKM) [19] was used to normalize the gene transcript levels.Transcript abundance of all genes in different samples was estimated by calculating the read density.DEGs were identi ed among fourteen libraries using EBSeq software [20].A threshold for false discovery rate (FDR) <0.01 and an absolute value of log2 fold-change (FC)≥1 were used to identify signi cant expression changes.

qRT-PCR analysis
In order to verify the expression abundance obtained by transcriptome sequencing, 15 genes were selected for qRT-PCR analysis in different samples.All the primers used for qRT-PCR were listed in Supplementary Table S1.cDNA synthesis and qRT-PCR analysis were performed as previously described [16].The relative expression level was calculated using 2 -ΔΔt method and standardized to the constitutive expression level of Actin [21].Three biological replicates were used for each gene.

Functional annotation
To acquire the most comprehensive annotation, functional annotations were performed by sequences similarity search against ve public databases.All genes were aligned with NCBI non-redundant protein database (Nr) [22], Gene ontology (GO) [23], Clusters of Orthologous Groups (COG) [24], SWISS-PROT [25] and the Kyoto Encyclopedia of Genes and Genomes databases (KEGG) [26].Blast2GO was used to obtain GO functional categories [27].

Selection of genes involved in speci c pathways and heatmap analysis
Glyphosate is proved to speci cally inhibit the synthesis of phenylalanine, tyrosine and tryptophan in the shikimate pathway, which serve as the precursors of numerous compounds such as pigments, alkaloids, hormones and cell wall components [40].Based on the annotation information in NCBI database, DEGs involved in some speci c pathways were selected for detailed analysis, and Heatmap analysis of DEGs was performed by BMKCloud (https://www.biocloud.net/).

Phenotypic and glyphosate tolerance analysis
To determine the different tolerance to glyphosate, seedlings of ZH10, MD12, ZH10-6 and HJ698 were treated with RoundupÒ at the dose of 0.18 g a.e./m 2 .There was no detectable growth difference at 12 hpt in these four genotypes after treatment.As the observation time increased (24, 72, 120 hpt), NT soybeans ZH10 and MD12 had visible foliar damage such as leaf roll, leaf chlorosis, and even leaf dehydration, whereas only leaf chlorosis could be observed in GM soybean HJ698 at 120 hpt but not in GM soybean ZH10-6.These results indicated that ZH10-6 and HJ698 showed tolerance to glyphosate and ZH10-6 displayed better growth condition after glyphosate treatment (Fig. 1).Interestingly, the "yellow ashing" phenomenon was only observed in terminal lea ets of HJ698 (Fig. 1).In other words, overexpression of EPSPS contributed to improve glyphosate tolerance in soybean, while overexpression of EPSPS and GAT conferred higher tolerance to glyphosate.

RNA-seq analysis
To compare the whole transcriptome response to glyphosate between these soybean genotypes, RNA samples were isolated from terminal lea ets after glyphosate treatment (0, 12 and 72 hpt), respectively.14 cDNA libraries were constructed for RNA-seq with two biological replications each for two time points (ZH10-6_12 hpt and HJ698_12 hpt).A total of 388 million raw reads were obtained and 90.72 Gb sequencing data was generated, with at least 4.8 Gb for each sample on average after data ltering.The average ratio of clean reads and GC content were about 91.5% and 49.8%, respectively.More than 92 % of the reads had Phred-like quality scores at Q30 level (Supplementary Table.S2), exhibiting that the sequencing data was high quality.To detect the distribution of generated sequence, all clean reads from each library were mapped to reference genome (Supplementary Fig. S1).After genome mapping with TopHat [28], 30.5 to 57.9 million clean reads was successfully mapped to the reference genome, and the average mapping ratio was about 72.7% (Supplementary Table.S2).Pearson's correlation coe cient (r 2 ) is an important index to re ect the relationships among the two biological replicates [29], and the value of r 2 of two replicated libraries was 0.9928 and 0.9836 in ZH10-6_12h and HJ698_12h, respectively, indicating that the RNA-seq data were reliable and reproducible.

Identi cation of DEGs
FPKM value of each gene was calculated for evaluating gene expression levels and all DEGs were identi ed (Supplementary Table.S3).To verify the reliability of RNA-Seq, 15 DEGs were randomly selected for qPCR assays.There was the same trend with less discrepancy between qPCR and RNA-seq results, indicating that the transcriptome data are reliable (Supplementary Table.S4).A total of 2,946 DEGs (1,365 up-regulated and 1,581 down-regulated genes) and 2,896 DEGs (1,119 up-regulated and 1,777 down-regulated genes) were identi ed in ZH10 and MD12 at 12 hpt, respectively (Fig. 3).In addition, 6,972 DEGs (3175 up-regulated and 3797 down-regulated genes) and 4,423 DEGs (1959 upregulated and 2464 down-regulated genes) in ZH10 and MD12 at 72 hpt respectively (Fig. 3).4,711 DEGs (1893 up-regulated and 2818 down-regulated genes) and 4,953 DEGs (2833 up-regulated and 2120 downregulated genes) were identi ed in HJ698 only expressed EPSPS at 12 and 72 hpt, respectively.7,820 DEGs (3319 up-regulated and 4501 down-regulated genes) and 2,289 DEGs (1103 up-regulated and 1186 down-regulated genes) were identi ed in ZH10-6 at 12 and 72 hpt, respectively (Fig. 3).Compared to 12 hpt, the numbers of DEGs slightly increased in HJ698 at 72 hpt while the numbers of DEGs signi cantly reduced in ZH10-6 at 72 hpt.

Similar response to glyphosate in two NT soybeans
Both GO and KEGG enrichment analysis were used for identi cation of key genes and pathways in response to glyphosate.GO enrichment analysis showed that a large number of DEGs were enriched in response to chitin (225 genes), response to wounding (195 genes), response to water deprivation (193 genes), and defense response to fungus (191 genes) terms in ZH10_12 hpt (Supplementary Fig. S2a), while a large number of DEGs were enriched in oxidation-reduction process (537 genes), chloroplast stroma (388 genes), chloroplast envelope (368 genes), chloroplast thylakoid membrane (292 genes) and chlorophyll biosynthetic process (115 genes) terms in ZH10_72 hpt (Supplementary Fig. S2b).Similar to ZH10, DEGs were mainly associated with plant stress and defense responses in MD12_12 hpt, while associated with photosynthesis and chlorophyll biosynthesis pathways in MD12_72 hpt (Supplementary Fig. S2c and S2d).KEGG analysis revealed that DEGs of ZH10_12 hpt were signi cantly grouped into plant hormone signal transduction, phenylpropanoid biosynthesis and plant-pathogen interaction categories.Among them, plant hormone signal transduction and plant-pathogen interaction pathways were associated with plant stress and defense response (Supplementary Fig. S4a).However, a different clustering tendency was identi ed at ZH10_72 hpt.DEGs grouped into 20 clusters and mainly contained porphyrin and chlorophyll metabolism, carotenoid biosynthesis, pentose phosphate pathway, glycolysis/gluconeogenesis, carbon xation in photosynthetic organisms, photosynthesis-antenna proteins, photosynthesis pathways (Supplementary Fig. S4b).These pathways were mainly related to photosynthesis, chlorophyll metabolism, and some pathways such as pentose phosphate pathway, glycolysis and phenylpropanoid biosynthesis were located in the upstream or downstream of shikimate acid pathway.Similarly, KEGG analysis revealed that DEGs of MD12_12 hpt were signi cantly associated with plant stress response and defense response pathways, and DEGs of MD12_72 hpt were mainly enriched in chlorophyll biosynthesis, photosynthesis, and shikimate-related pathways (Supplementary Fig. S4c and S4d).Moreover, a few DEGs of MD12_12 hpt were found to involve in photosynthesis and carotenoid biosynthesis pathways (Supplementary Fig. S4c).
Different response to glyphosate in two GM soybeans GO enrichment analysis revealed that DEGs of HJ698_12 hpt were classi ed into 20 clusters and mainly involved in plant stress response and plant defense response pathways, including oxidation-reduction process (323 genes), response to water deprivation (269 genes), response to chitin (258 genes), response to wounding (257 genes), and defense response to fungus (257 genes) (Supplementary Fig. S3a).GO enrichment analysis revealed that DEGs of HJ698_72 hpt were also mainly grouped into plant stress response and defense response pathways, and only a few DEGs were involved in chlorophyll photosynthesis and shikimate-related pathway (Supplementary Fig. S3b).Unlike two NT soybeans and HJ698, a signi cant difference was found in ZH10-6_12 hpt.GO enrichment analysis revealed that DEGs of ZH10-6_12 hpt were mainly involved in function of chlorophyll and photosynthesis pathways, including chloroplast stroma (369 genes), chloroplast envelope (340 genes), and chloroplast thylakoid membrane (206 genes) (Supplementary Fig. S3c).Additionally, DEGs of ZH10-6_72 hpt were mainly enriched into plants response to stress and defense pathways, which contained response to chitin (142 genes), response to water deprivation (138 genes), response to wounding (136 genes), defense response (133 genes), response to cold (128 genes) (Supplementary Fig. S3d).KEGG analysis revealed that DEGs were clustered into similar function categories both in HJ698 at 12 hpt and 72 hpt.These genes were mainly involved in carbon metabolism, biosynthesis of amino acids, plantpathogen interaction, avonoid biosynthesis, phenylpropanoid biosynthesis, pentose phosphate pathway, glycolysis/gluconeogenesis, and photosynthesis (Supplementary Fig. S5a and 5b).KEGG analysis revealed that DEGs of ZH10-6_12 hpt were mainly involved in photosynthesis and energy metabolism pathway, including carbon metabolism, carbon xation in photosynthesis organisms, photosynthesis, TCA cycle, glycolysis/gluconeogenesis, and pentose phosphate pathway (Supplementary Fig. S5c), while DEGs of ZH10-6_72 hpt were signi cantly associated with plant hormone signal transduction, phenylalanine metabolism, circadian rhythm-plant, fatty acid elongation and avonoid biosynthesis, but chlorophyll biosynthesis, photosynthesis and shikimate pathway was not found in ZH10-6_72 hpt (Supplementary Fig. S5d).

DEGs involved in shikimate acid biosynthetic pathway
Since glyphosate is proved to interfere with the shikimate acid biosynthetic pathway, genes involved in this pathway were selected for further analysise.18 and 8 DEGs involved in shikimate acid pathway were identi ed at 12 hpt in ZH10 and MD12 respectively, but this number was obviously increased into 33 and 18 at 72 hpt.Compared to NT soybeans, a relative gentle increasing trend of DEGs was observed in HJ698 (20 and 22 DEGs were identi ed at 12 and 72 hpt, respectively).However, a completely opposite change trend of DEGs was observed in ZH10-6, in which 32 DEGs involved in shikimate acid pathway were identi ed at 12 hpt while only 6 were identi ed at 72 hpt.Meanwhile, to compare the changes of shikimate acid biosynthetic pathway among these soybeans, a heat map was generated to present the transcript abundance for all DEGs.Heatmap culster analysis revealed that expression level of many genes were continuously increased or decreased in NT soybeans and slightly altered in HJ698 as the extension of treat time.As to genes in ZH10-6, the trend of gene expression was totally different.The expression level of these genes was only signi cantly altered at 12 hpt, but it tended to recover to levels of 0 hpt at 72 hpt (Fig. 4).These results indicated that shikimate acid biosynthetic pathway was associated with different phenotypic variation among these soybeans.

DEGs involved in herbicide targeted cross-pathways
Interestingly, inhibition of shikimate pathway by glyphosate can speci cally in uence the expression of genes not only in target pathway but also cross-pathways, and the action modes of herbicides mainly include photosystem, chlorophyll biosynthesis, synthetic auxin, microtubule and cellulose biosynthesis pathways [18].Signi cantly, the "yellow ashing" phenomenon was also observed in glyphosate-tolerant soybean HJ698 after glyphosate application (Fig. 1)., Previous study revealed that the decrease of chlorophyll content is the main reason for this symptom [16].Detailed analysis of genes related to photosystem and chlorophyll biosynthesis pathways showed that the expression level of many DEGs at 72 hpt was obviously altered than that at 12 hpt in NT soybeans, and slightly altered in HJ698.However, the expression level of these genes was signi cantly altered at 12hpt but it tended to recover to levels of 0 hpt at 72 hpt in ZH10-6(Fig.5).
DEGs involved in synthetic auxin, microtubule and cellulose biosynthesis pathways.A total of 117 DEGs related to synthetic auxin and 43 DEGs involved in microtubule and cellulose biosynthesis response to glyphosate in four soybean genotypes were identi ed.The number of DEGs related to synthetic auxin, microtubule and cellulose biosynthesis at 72 hpt was obviously increased in HJ698 and NT soybeans compared to 12 hpt.However a reverse change trends of No. of DEGs related to synthetic auxin and microtubule and cellulose biosynthesis was observed in ZH10-6, 60 and 21 DEGs were identi ed at 12 hpt while only 12 and 6 DEGs were found at 72 hpt.Interestingly, the variation tendency of expression abundances of DEGs was obviously strengthened in NT soybeans and HJ698 as the increasing of treat time.But the abundance level at 72 hpt was inclined to return to the level at 0 hpt (Supplementary Fig. S6, S7 and S8).This result was consistent with the analysis of expression abundance of DEGs related to shikimate pathway, photosystem and chlorophyll biosynthesis pathways.These results showed that photosystem, chlorophyll biosynthesis, synthetic auxin, microtubule and cellulose biosynthesis were also associated with phenotypic variation among these soybeans after glyphosate application.

Discussion
In this study, the transcriptome analysis of four soybean genotypes in response to glyphosate were revealed by RNA-seq, and these transcriptional changes will contribute to the comprehensive understanding of DEGs or pathways involved in response to glyphosate.Similar to pathogen attacks, plant defense and stress response mechanism will be activated when glyphosate exposure, and oxidation-reduction process was indispensable defense strategy.GO enrichment analysis showed that DEGs involved in plants stress response and plant defense were mainly in NT soybeans at 12 hpt, in HJ698 at 12 and 72 hpt, and in ZH10-6 at 12 hpt (Supplementary Fig. S2 and S3).Surprisingly, DEGs associated with plant stress response and plant defense were not enriched in ZH10-6 at 12 hpt (Supplementary Fig S3c).
Due to the non-selective characteristic of glyphosate, it is easily accumulated in plants through glyphosate spray.The uptake of glyphosate is a biphasic process which contains glyphosate penetration through the cuticle and glyphosate translocation into the symplast with phosphate transporters [30,31].
Excessive glyphosate accumulation led to increased biological damage effects, which can be detected in NT soybeans and HJ698 (Fig 1).Visible foliar damages such as"yellow ashing" phenomenon, leaf roll and even leaf dehydration were also observed in two NT soybeans ZH10 and MD12 (Fig 1).Some intriguing studies showed that the "yellow ashing" symptom was observed in non-transgenic soybeans and GM soybean only expressed with EPSPS gene after glyphosate application [7,9,16].It is worth noting that damage will easily take place in the actively growing tissues of plants if glyphosate is not su ciently degraded.The applied glyphosate was accumulated in the meristematic region of TN and GM cottons, and overexpression of EPSPS and GAT genes conferred lower glyphosate residues by glyphosate acetylation in cotton [8].Similar to the reported results, the "yellow ashing" phenomenon was observed in the terminal lea ets of HJ698, while not in the terminal lea ets of ZH10-6 when introduction of GAT gene (Fig 1).Therefore, the residual dose of glyphosate after glyphosate acetylation was not enough to cause a toxic effect in ZH10-6.Importantly, the decrease of chlorophyll content is closely related to the "yellow ashing" phenomenon, then affects the photosynthetic process in plant [10,11].However, DEGs associated with function of chlorophyll and photosynthesis pathways were mainly enriched in ZH10-6 at 12 hpt after glyphosate application.Compared to 12 hpt, the numbers of DEGs related to photosystem and chlorophyll biosynthesis pathways were obviously increased in two NT soybeans, while slightly increased in HJ698, but signi cantly decreased in ZH10-6 at 72 hpt (Supplementary Fig. S2 and S3).
Inhibition of plant growth is another typical symptom after glyphosate application [35].Additionally, synthetic auxin, microtubule and cellulose biosynthesis are critical for plant growth and development process [18].As was shown in Fig. 1, the plant height was signi cantly inhibited in two NT soybeans at 120 hpt.Meanwhile, compared to 12 hpt, DEGs involved in synthetic auxin and microtubule and cellulose biosynthesis were obviously increased in HJ698 and NT soybeans at 72 hpt, but signi cantly decreased in ZH10-6 at 72 hpt.(Supplementary Fig. S6 and S7).These results showed that overexpression of EPSPS gene or/and GAT gene improved the genes expression related to synthetic auxin, and microtubule and cellulose biosynthesis pathways for keeping plant growth.Phenylpropanoid pathway plays an important role in plant growth and development as well as in biotic and abiotic stress responses, and oxidative stress responses [17,36,37,38], and phenylpropanoid pathway also can be found in ZH10, MD12, ZH10 and MD12 (Supplementary Fig. S4 and S5).These results showed that phenylpropanoid pathway, synthetic auxin, microtubule and cellulose biosynthesis were also associated with phenotypic variation among these soybeans.Additionally, the numbers of DEGs involved in shikimate pathway were increased, and their contents of shikimate acid also accumulated rapidly in NT plant [33,34].Compared to HJ698 and NT soybeans, the number and expression level of DEGs related to shikimate pathway was obviously reduced in ZH10-6 as the extension of glyphosate post-treatment (Fig. 4).Previous study showed that numbers of DEGs involved in shikimate pathway was not found in NT soybean at 6 hpt after glyphosate treatment, but subsequently increased at 24 hpt and 72 hpt [18].Maybe it would take some time for glyphosate uptake before inhibition of shikimate pathway.Surprisingly, the expression trend of 10 DEGs involved in shikimate pathway was completely opposite between ZH10-6 and HJ698 at 12 hpt (Supplementary Fig. S8).These results indicated that the GAT gene might work at the early stage under glyphosate treatment, and subsequently EPSPS gene will be activated when glyphosate uptake.
Taken together, introduction of EPSPS and GAT genes conferred higher tolerance to glyphosate and showed better growth condition, and its comparative transcriptome changes were detected after glyphosate treatment, which will be valuable for further efforts to elucidate molecular mechanisms underlying glyphosate tolerance or glyphosate detoxication.

Abbreviations
Clusters of Orthologous Groups; DEGs: Differentially expressed genes; EPSPS: 5enolypyruvylshikimate-3-phosphate synthase; FDR: false discovery rate; FPKM: Fragments per Kilobase of transcript per Million fragments mapped; GAT: Glyphosate N-acetyltransferase; GM: Genetically modi ed; GO: Gene ontology; GOX: Glyphosate oxidoreductase; KEGG: Kyoto Encyclopedia of Genes and Genomes databases; NT: No-transgenic; RNA-Seq: Transcriptome sequencing.Declarations Ethical approval and consent to participateSoybeans used in this study were permitted to do scienti c research.The experimental research on soybean was complied with Chinese legislation and eld studies were in accordance with guidelines of Institute of Crop Science, Chinese Academy of Agricultural Sciences.All the experiments were not involved in any endangered or protected species.

Figures
Figures