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. /m2. 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 flashing” phenomenon was only observed in terminal leaflets 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 leaflets 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 filtering. 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 coefficient (r2) is an important index to reflect the relationships among the two biological replicates [29], and the value of r2 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.
Identification of DEGs
FPKM value of each gene was calculated for evaluating gene expression levels and all DEGs were identified (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 identified 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 up-regulated 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 down-regulated genes) were identified 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 identified 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 significantly reduced in ZH10-6 at 72 hpt.
Similar response to glyphosate in two NT soybeans
Both GO and KEGG enrichment analysis were used for identification 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 significantly 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 identified at ZH10_72 hpt. DEGs grouped into 20 clusters and mainly contained porphyrin and chlorophyll metabolism, carotenoid biosynthesis, pentose phosphate pathway, glycolysis/gluconeogenesis, carbon fixation 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 significantly 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 classified 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 significant 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, plant-pathogen interaction, flavonoid 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 fixation in photosynthesis organisms, photosynthesis, TCA cycle, glycolysis/gluconeogenesis, and pentose phosphate pathway (Supplementary Fig. S5c), while DEGs of ZH10-6_72 hpt were significantly associated with plant hormone signal transduction, phenylalanine metabolism, circadian rhythm-plant, fatty acid elongation and flavonoid 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 identified 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 identified 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 identified at 12 hpt while only 6 were identified 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 significantly 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 specifically influence 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]. Significantly, the “yellow flashing” 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 significantly 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 identified. 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 identified 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.