ABSTRACT
ABSTRACT
Genetically modified (GM) organisms have been developed for decades. However, unintended effects are the main concerns of safety assessment that needs to be carefully investigated. Here, eight varieties of GM rice that were developed in China were selected to assess the unintended effects through transcriptome and metabolism. There are 2892–8758 differentially expressed genes (DEGs) and 7–50 metabolites at significant level between GM varieties and their isogenic counterparts, which were far fewer than that between traditional rice varieties. The function enrichment analysis showed altered transcription in stress-related pathway and starch and sucrose metabolism. DEGs shared among eight GM samples constitute less than 1% of the genes in the genome, and none of them is reported more than four times. The insertion effect on the nearby gene expression and the associated metabolism is only restricted to 50 genes. All the results provide a comprehensive analysis of unintended effects and indication of difference in Chinese transgenic rice based on their backgrounds, transformation, and insertion elements.
INTRODUCTION
Genetically modified (GM) crops have been developed rapidly since the mid-eighties, possessing 185.1 million hectares in. 2016.12 The top ten countries grow more than 1 million hectares each.1 GM crops are manipulated by inserting target DNA fragments into the host genome to improve traits, which has been developing into insect-resistant, herbicide-resistant, quality improvement, and stacked traits. For this preference, GM crops have become the possible alternative to improve grain yield and reduce labor cost.3
The public and regulators’ concerns on GM crops mainly focus on whether GM crops are subjected to strict scrutiny for their safety.4 It has been agreed internationally through countries and commissions to ensure an substantial equivalence between GM crops and traditional varieties, which is approved internationally.5–8 The USA has published its research on effects of GM crops on human health, environment, agronomy, and economy; they found no adverse effects attributed to genetic engineering in human population.9 However, the majority of commercialized GM crops are created by Agrobacterium-mediated gene transfer10–12 and particle bombardment13–15 that result in random insertion4 and new proteins, new metabolites or altered levels of existing metabolites that could compromise safety.8 It needs to be a comprehensive analyzing method to compare the GM crops or products to their counterparts from parental or near isogenic lines.16,17
The role of emerging “omics” technologies in the assessment of unintended effects was recently proposed and has been applied for the reconstruction of genome-scale networks for model organisms in the commercialized GM crops,16,18,19 such as maize, soybean, canola, potato, and stacked GM maize variety.20–26 Unfortunately, current studies are mostly based on individual events, and it is difficult to evaluate whether there is uniform change in genome-scale.
Rice (Oryza sativa L.) is one of the most important food crops, providing food for over half the world’s population.27 Researches on transgenic rice have been extensively carried out, but few are focused on the safety management.28–30 Although seven transgenic rice cultivars are reported to be approved for cultivation or food and feed use in 2017,1 most of them are not commercialized across the world. Rice is one of the major crops in China, and transgenic insect-resistant rice TT51-1 (synonym BT63) that harbors a hybrid Cry1Ab/Ac gene was granted a safety certificate for commercial use in 2009 as the first transgenic food crop.31–33 However, limited published results guided us for better understanding the unintended effects of GM rice.30,34,35
In this study, eight genetically engineered rice lines in China and their corresponding wild types were collected. Changes of GM plants at whole genome level and the insertion site for transcriptions and metabolites compared to their counterparts were explored. We evaluated the changes in gene expression and metabolites between GM rice and different traditional maize varieties and between traditional maize varieties. This work provided a comprehensive analysis to evaluate the unintended effects of GM rice and a safety analyzing result of GM rice developed in China.
METHODS
Plant Materials
The experimental transgenic materials and their corresponding counterparts are described in Table 1, which is proved to be genetically stable by providers. GM1 and wild type 1 (WT1) were provided by Prof. Daichang Yang (Wuhan University). GM3, GM4, and WT3 were provided by Prof. Huachen Yu Chinese Academy of Sciences. GM5/WT4 and GM6/WT4 were provided by Prof. Qifa Zhang (Huazhong Agriculture University). GM8 were provided by Prof. Wensheng Huang (Chinese Academy of Inspection and Quarantine). GM2, GM7, WT2, and WT5 were stored at our laboratory. GM1 and GM2 were special because they both possess unique gene inserted to unique parental line, which is completely different from other materials.
TABLE 1. Sample information included in this study.
Samples Preparation
Seeds of rice were surface sterilized with 10% H2O2 (v/v) for 10 min, rinsed thoroughly with distilled water, and germinated on moist filter paper for 3 days in an incubator at 25°C. After germination, the seeds were then transferred to a net floating on 0.5 mM CaCl2 in a plastic container. After 6 days, the seedlings were transferred to a 10-L plastic pot in a growth chamber at 400 μM photons/m2 s and 28°C/18°C under a 14/10-h light/dark regime. The nutrient solution for rice was half-strength Kimura B solution containing the macronutrients, including 0.18 mM (NH4)2SO4, 0.27 mM MgSO4·7H2O, 0.09 mM KNO3, 0.18 mM Ca(NO3)2·4H2O, and 0.09 mM KH2PO4 and the micronutrients, including 20 μM NaEDTAFe·3H2O, 6.7 μM MnCl2·4H2O, 9.4 μM H3BO3, 0.015 μM (NH4)6Mo7O24·4H2O, 0.15 μM ZnSO4·7H2O, and 0.16 mM CuSO4·5H2O. The nutrient solution was prepared with distilled water, aerated daily, and renewed every 3 days. The leaves were harvested at the trefoil stage, frozen with liquid nitrogen for RNA extraction. Each sample was prepared in triplicate for sampling pooling and sequencing.
RNA-Seq and Bioinformatic Analysis
Total RNA was extracted using a protocol based on the Trizol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) and purified with the Qiagen RNeasy MiniElute Cleanup kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Extracted RNA was then quantified by a NanoDrop N2000 spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA). OD 260/280 nm absorption ratios (mean = 2.02 ± 0.02) were used to confirm the integrity and purity of the RNA samples. The library for sequencing was prepared using the Illumina mRNA-Seq 8-Sample Prep Kit following the manufacturer’s instructions (Illumina, USA). The 125 base paired-end sequencing was performed on the Illumina Hiseq-2000 platform (Illumina, USA). The sequenced reads were further quality-filtered using FastQC 0.11.5 tools to remove poor-quality sequences. Clean reads were mapped to the Oryza Sativa reference genome IRGSP-1.0.22 (ftp://ftp.ensemblgenomes.org/pub/release-38/plants/fasta/oryza_sativa/) using TopHat 2.1.0 with the default parameters. The normalized expression level (fragments per kilobase million) for each gene was then estimated using Cufflinks 2.2.0 based on the alignment from the last step.42 We used analysis of variance to detect the DEGs for each pair of transgenic rice line and their isogenic wild type line through in-house SAS scripts. DEGs among isogenic wild type were also calculated through this method. Principle component analysis (PCA), hierarchical clustering, and K-means analysis to compare the difference between GM and WT lines were performed through in-house SAS scripts.
Metabolite Profiling
The metabolomic profiling analysis was performed by Metabolon (Durham, NC, USA).43 The metabolomic platforms consisted of three independent platforms: ultrahigh performance liquid chromatography/tandem-mass spectrometry (UPLC/MS/MS2) optimized for basic species, UPLC/MS/MS2 optimized for acidic species, and gas chromatography/mass spectrometry (GC/MS) for volatile species. Three biological replicates were provided for each sample.
Pathway Mapping and Analysis
Genes with Ensembl ID were integrated and mapped on gene ontology (GO) terms (http://systemsbiology.cau.edu.cn/agriGOv2/download/871_slimGO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (http://structuralbiology.cau.edu.cn/PlantGSEA/database/Osa_KEGG) through in-house SAS scripts. The Pathview package was used to perform the KEGG pathway analysis of metabolites.44 Metabolites of each sample at significant level were integrated and linked with related protein and genes through STITCH website 5.0 (http://stitch.embl.de/cgi/).
Insertion Analysis
Flanking sequence of each GM line was first obtained from publications and Chinese applications (Table 1). Their locations were determined by blasting analysis of flanking sequence to Oryza Sativa reference genome. The coding genes surrounding the insertion site were retrieved and grouped according to the distance to the insertion site. Each cluster contained several numbers of genes. The average P value of genes in one cluster was calculated. These five P values were then used to score the inversion number to demonstrate the insertion influence on nearby genes. Such analysis was performed five times, as numbers of coding genes contained in one cluster were 1, 3, 5, 10, 20, 30, and 50, respectively. A low inversion number, which should be derived from the consecutive increase of average P value in the further adjacent clusters, would indicate the existence of hierarchical continuous effect due to gene insertion. Negative control was calculated based on the corresponding insertion site and the nearby transcriptome of non-GM materials.
qRT-PCR
Total RNAs were isolated with a Trizol reagent (Invitrogen, CA, USA). First-strand cDNA was synthesized using a ReverTra Ace® qPCR RT Kit (TOYOBO, Japan), and PCR reaction was performed with a SYBR Green® Realtime PCR Master Mix (TOYOBO, Japan) on a LightCycler 480 II (Roche, CH). Each reaction mix contains 2 μL cDNA, 0.3 μM both primers, 10 μL Master Mix, and ddH2O was added to form a total volume of 20 μL. The PCR reaction condition was as follows: 95°C for 1 min, following 95°C for 10 s, 60°C for 60 s for 40 cycles. A fluorescence signal was detected at the end of each cycle. Three replicates were arranged for every data point. Actin was used as an internal control.
RESULTS
GM Rice Lines Show Less Differences to Wild Counterparts than WT Comparisons
We obtained more than 903 million 125 bp pair-end reads that passed the quality filters in these 13 samples; each of which is analyzed by triplicate sequencing data. The metabolism data obtained from GC–MS and LC–MS are also satisfied for further analysis (5: Figure S1 and 6: Figure S2).
A total of 218 metabolites and 35,678 genes were quantified across 13 lines, and the DEGs have been selectively verified through q-PCR experiment (Table S4). There are 2892–8758 DEGs (P < 0.05) between GM variety and their isogenic counterparts (Table 2), which were far fewer than that between traditional rice varieties (t-test, DF = 16, t = 3.31, P = 0.004). The same trend was also found at the level of metabolites (Table 2), which means that only 7–50 metabolites at significant level were found in WG (differential analysis between genetically modified rices and their corresponding counterparts) comparisons (t-test, DF = 16, t = 2.77, P = 0.014). After subtracting those overlapped genes that existed in both WW (differential analysis between two isogenic counterparts in this research) and WG comparisons, the number of DEGs (unique mRNAs in Table 2) was decreased to 164–667. These results indicate that transcriptional turbulences caused by transgenic constructs are much less than that by genetic and environmental influence. This is consistent with the previous findings that environment and/or naturally species variations affect gene expression more strongly than the genetic modification.21,45
TABLE 2. Significant gene expressions and metabolites summary of all GM rice varieties and their corresponding counterparts.
PCA was further used to show the differences between the WG and the WW comparisons (Fig. 1). About 45.6% of the total variance was explained by the first three PCs (PC1, 26%; PC2, 13%; and PC3, 5%), suggesting that these wild types are moderately separated by gene expression patterns (Fig. 1A). GM samples tend to be grouped with the respective non-GM counterparts in plots (Fig. 1B). Analysis on metabolites shows similar trends as the expression data (Fig. 1C,D). Hierarchical clustering showed that each GM sample and its corresponding counterpart belong to the same class (Fig. 1). This is the first comprehensive analysis on major Chinese GM rice events that GM events seem to introduce less variation on the transcripts and metabolites than that between non-GM counterparts.
Published online:
08 April 2019Figure 1. PCA and clustering analysis of transcripts (A, B, and E) and metabolites (C, D, and F) of GM and WT samples, respectively.
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/kgmc20/2019/kgmc20.v010.i01/21645698.2019.1598215/20190610/images/medium/kgmc_a_1598215_f0001_oc.jpg"})
Figure 1. PCA and clustering analysis of transcripts (A, B, and E) and metabolites (C, D, and F) of GM and WT samples, respectively.
Few Overlapped DEGs, Metabolites, and Enriched Functions were Found Between GM Lines
For unique DEGs reported by each WG comparison, less than 15% were reported two times. Sixty six and eight genes were reported three or four times, respectively (Table S1). The function of these genes mainly attributed to response to external stimulus or stress and transformation is obviously one of these stresses (Table 3). The metabolites analysis offered the similar conclusions, which means that nine overlapped metabolites that reported three times at significant concentration only influenced amino acid metabolism between WG comparisons (Table S1).
TABLE 4. K-means analysis of WG comparisons.
TABLE 3. List of gene IDs that reported four times and KEGG compound IDs that reported three times.
The K-means analysis further revealed that WG groups are close to each other based on their parental genotype and the insertion element (Table 4), such as WG7 and WG8 that belong to the same class because of their similar parent origin. However, GM5 and GM6 are also based on the same parent background but belong to two different groups, indicating that effects caused by transgene are not solely parent dependent. Moreover, GM5 and GM6 are based on the same isogenic wild type counterpart (Minghui 63) and shared 18 DEGs and no metabolites at significant level, which is mainly related to RNA splicing and biosynthesis of secondary metabolites (Table S2). GM3 and GM4 were developed using the same construct and transformation method (pFHBT1) and shared four DEGs and no metabolites. Difference between GM rice and non-GM counterparts is related to insertion element, background, and other factors and should be evaluated case-by-case.
The GO analysis mainly contributed to the process of macromolecule metabolism, DNA and RNA biosynthesis, and plant stress/defense signaling. The representative GO and KEGG enrichment with a cutoff P < 0.05 is shown in Table 5 and Table S3. Most GO overlaps contributed to the biological and energy metabolism, which appeared five times. GM3 was distantly placed with WT3 in PCA analysis (Fig. 1); however, the GO analysis of WG3 was mainly enriched with defense response. Of all the six KEGG pathways, metabolic pathways exist in all eight GM events. Three GM events (WG2, WG5, and WG8) showed transcriptional modification in starch and sucrose metabolism that is related to plant defense response to the biotic or environment stress.46–48 These results tend to indicate that the acquisition of desired traits is accompanied by transcript levels of regulation of transcriptional process that is associated with environmental stress responses.49–52
TABLE 5. GO and KEGG enrichment of significantly induced and repressed genes for all WG comparisons.
Metabolites are Not Linked with DEGs in Response to Transgene
In order to investigate the effect of transgene on the entire metabolic pathways that functioned from transcript to metabolites, we linked the metabolites in each WG comparison to the corresponding coding genes and evaluated their expression. Each WG comparison has no more than two DEGs that are directly linked to metabolism profiling; most of which were related to citrate cycle and amino acid metabolism (Table 6). Additionally, WG1, WG3, and WG6 all showed variation on the defense response, mainly focused on the transmembrane transportation, metal tolerance, and transcription activation.
TABLE 6. Differentially expressed genes that associated with metabolites and their functions.
Insertion Analysis Suggest only Short Range Affections
By blasting analysis of the recorded flanking sequences from each transgenic event to Oryza Sativa reference genome, we identified the insertion sites of each transgenic event in the gnome. Results showed that the insertion sites of WG1, WG5, WG7, and WG8 were located at the coding sequence, while all others are located on the intergenic region. The insertion site was also analyzed to determine if it locates in the high-expression region or gene enrichment region, as these regions seemed much significant on the gene expression and the corresponding metabolic regulation. By randomly selecting sampling gene expression for 5000 times, the insertion site of GM5 was determined to be located on the relatively high-expression region.
To test whether insertion results in a server disturb on gene expression for the nearby genes than the remote genes, we calculated the inversion number (the cardinality of inversion set of a set of numbers) for each sample (step size in Table 7). For each test in the WG comparison, a higher P value indicates less changes at transcriptional level. To describe the trends of ripple effect, we calculated the inversion number based on P values for the adjacent clusters. We calculated the inversion number of each step size. As shown in Table 6, we did not observe zero for any step size, and we did not observe a clear trend of increasing inversion number with increasing step size either. It suggests that the ripple effect might not exist. The insertion event would only bring short range effect. By using random WW comparison as negative control, we recalculated the inversion number and found the similar result with that from the canonical WG comparison, which confirms the sporadic effect on the nearby genes. We next test the short-range effect by checking if the P value of the first cluster (P1) is the smallest in that of all the five clusters. By using the results from 100 permutations as negative control, we found GM3 and GM4, but not GM5-exhibited influence on up to three nearest genes, while all other insertion events failed to pass the threshold from the null distribution of inversion number. Altogether, the insertion events did not impose long-range effects on the gene expression; even if it located on the high expression region. Even in few cases, short-range effects are observed, and it seems not to be a common, but case-specific, pattern.
TABLE 7. Inversion number calculation of each stepsize between GM materials and negative control.
DISCUSSION
Unintended effects of GM crops are the main concern of administrators and citizens after their intended effects have been thoroughly investigated. Many efforts have been made to test if GM crops are substantially equivalent as safe non-GM comparators, following a large body of high-quality compositional data (around 50–150 analytes) that need to be determined according to principles outlined in the Organization for Economic Cooperation and Development.7 This concern promoted the use of modern technologies, especially an unbiased high-throughput “omics” technology to validate the substantial equivalence.16 In these transgenic rice varieties, Bt genes and the expressed proteins are not native to plants and exert no known metabolic activity in plant, and so is the OsrHSA recombinant protein. However, we still found eight DEGs that reported four times in WG comparisons, especially in insect-resistant Bt rice (Table 3). Of these unintended DEGs, two are related to starch and sucrose metabolism, representing sucrose-phosphatase (OS05G0144900) and trehalose-phosphate phosphatase (OS06G0222100), respectively. These enzymes are mainly involved in the sugar metabolism and photosynthesis.53,54 Under drought stress conditions, the critical roles of sucrose, glycine-betaine, and trehalose in the starch and sucrose metabolism (ko00500) have been actively researched to understand the tolerance of plants to biotics, dehydration, and insect invasion.46–48 Other genes are related to plant response to abiotic response, such as OPR (OS08G0459600) that is functional in the metabolism of jasmonic acid biosynthesis,55 TFIIF (OS10G0188100) that is involved in the process of transcription,56 and ZIF (OS11G0135000) that involved in the signaling of plant immunity.57 Additionally, DEGs that were linked directly to metabolism profiling were mainly involved in amino acid metabolism, which could partly because of the increased demand of inserted protein biosynthesis.58 These DEGs expressed in GM rice should be further studied for their edible and environmental safety.
The environment was even shown to play a stronger effect in the protein, gene expression, and metabolite levels of the GM samples than gene modification by previous publications, which showed similar results for transcriptomics in maize, wheat, soybean, and potato.20–22,25,45,59 Our research on GM rice developed in China indicated less variation in gene expression and metabolites than non-GM varieties (Table 2). Only half of GM events shared DEGs without function enrichment so that gene insertion did not pose metabolic change on the rice physiology. WG2 and WG3 are associated with some stress-related mechanisms such as programmed cell death and apoptotic process.60–62 However, Bt rice is resistant to insects through making them lethal after taking Bt protein, which inevitably stimulate some defense pathways according to the previous report.63 High-level expression of Bt genes as well as maker genes and artificial sequence of gene regulators in this rice may act as a physiological sensor to disturb stress-related pathways. This observation has been proved in the assessment of GM2 and GM3 previously.58,63,64 These could possibly be common effects caused by Bt transgene. However, different DEGs and metabolites in significance existed in WG comparisons suggest that this effect may also be related to construct, background, and transformation method. Moreover, gene insertion-induced pathological variation occurred in WG3, but not in WG4 that was cultivated with GM3 and Zhenshan 97, providing an implication that crossover between GM and dominant rice that compromise the unintended effects to stabilize the genetics of inserted gene. Gene insertion is indicative of the regulation of stress-related genes and pathways that may act a role to protect plant from environmental changes and should fully be evaluated for plant genetic transformation technologies. This observation was also validated through linkage analysis between DEGs and metabolites at significant level, where nine candidate DEGs that linked metabolites at significant level attributes to metabolic process and stress-related pathways. This result also suggested that GM rice by inserting foreign genes was much like the process of exposure to abiotic stress, which induced the expression of defense genes for adaptation.46,47
This study provides an unbiased, comprehensive profiling of unintended effects in GM rice. Considering their breeding development, the extensive vetting is involved in the generation and selection of one or a few “elite” events, which artificially prevents the “bad” GM crops that might produce “bad” metabolites or unintended effects for food and feed use.26 On the other hand, it is also revealed in this study that difference of metabolite levels is less than that of transcripts level, no matter in GM or WT lines (Figure S3). This is also corroborated by previous studies on variation in human chromatin and primate transcript state.65,66 GM rice is genetically stable through molecular breeding and undergoing assessment to be approved for commercialization, and their unintended effects based on omics technologies should be evaluated further.
CONCLUSIONS
Transcriptomic and metabolic analysis of Chinese transgenic rice that are potential for further commercialization and their isogenic counterparts were conducted in this study. Changes of transgene insertion in rice exposed less difference than other rice varieties. Difference between GM rice and non-GM counterparts is limited to defense-related pathways and is related to insertion element, background, and transformation process. Gene insertion in the genome does not impose long-range effects on the gene expression, even if it located on the high expression region.
| Genetically modified (GM) materials | Event name | Wild type (WT) | Insert gene | Promoter | Terminator | Insertion site | References |
|---|---|---|---|---|---|---|---|
| GM1 | 114-7-2 | WT1 (Taipei 309) | Oryza sativa recombinant human serum albumin (OsrHSA) | Rice Gt13a promoter | NOS2 | Os04g0604600 | [36] |
| GM2 | KMD | WT2 (Xianyou 10) | Bt CrylAb(Hpt)1 | CaMV 35S | NOS | OS02G0620533-OS02G0620600 | [37] |
| GM3 | Huahui No.1 | WT3 (Minghui 63-Beijing) | Bt CrylAb/crylAc | Actin | NOS | OS10G0173800-OS10G0174548 | [38] |
| GM4 | Bt 63 (TT51-1) | WT3 (Minghui 63-Beijing) | Bt CrylAb/crylAc | Actin | NOS | OS10G0173800-OS10G0174548 | [33] |
| GM5 | T2A-1 | WT4 (Minghui 63-Wuhan) | Bt Cry2A | Maize ubiquitin promoter | NOS | Os12g0631600 | [39] |
| GM6 | T1C-19 | WT4 (Minghui 63-Wuhan) | Bt CrylC | Maize ubiquitin promoter | NOS | OS11G0124300-OS11G0124500 | [40] |
| GM7 | Kefeng-6 | WT5 (Minghui 86) | Bt CrylAc/CptII(Hpt)3 | Maize ubiquitin promoter | NOS | Os04g0648750 | [41] |
| GM8 | Kefeng-8 | WT5 (Minghui 86) | Bt CrylAc/CptII4 | Maize ubiquitin promoter | NOS | Os11g0430000 | [41] |
1Bt gene was from Bacillus thuringiensis to confer resistance to insects, such as Bt Cry1Ab, Bt Cry1Ac, Bt Cry2A, Bt Cry1Ab/1Ac, etc.
2NOS is the terminator of nopaline synthase gene.
3Hpt gene is hygromycin phosphotransferase gene from Streptomyces hygroscopicus.
4CptII is a cowpea trypsin inhibitor.
| WW comparisons1 | Metabolites | mRNAs | Ratio of metabolites | Ratio of mRNAs | WG comparisons2 | Metabolites | mRNAs | Ratio of metabolites | Ratio of mRNAs | Unique mRNAs |
|---|---|---|---|---|---|---|---|---|---|---|
| W12 | 74 | 9538 | 33.94% | 26.73% | WT1/GM1(WG1) | 33 | 2566 | 15.14% | 7.19% | 216 |
| W13 | 88 | 11454 | 40.37% | 32.10% | WT2/GM2(WG2) | 50 | 5174 | 22.94% | 14.50% | 298 |
| W14 | 54 | 12403 | 24.77% | 34.76% | WT3/GM3(WG3) | 15 | 2892 | 7.34% | 8.11% | 185 |
| W15 | 58 | 16390 | 26.61% | 45.94% | WT3/GM4(WG4) | 10 | 4977 | 4.59% | 13.95% | 667 |
| W23 | 56 | 3954 | 25.69% | 11.08% | WT4/GM5(WG5) | 7 | 8758 | 3.21% | 24.55% | 164 |
| W24 | 25 | 8413 | 11.47% | 23.58% | WT4/GM6(WG6) | 8 | 3681 | 3.67% | 10.32% | 167 |
| W25 | 19 | 9880 | 8.72% | 27.69% | WT5/GM7(WG7) | 23 | 7301 | 10.55% | 20.46% | 250 |
| W34 | 30 | 4072 | 13.76% | 11.41% | WT5/GM8(WG8) | 19 | 3027 | 8.72% | 8.48% | 396 |
| W35 | 52 | 9651 | 23.85% | 27.05% | ||||||
| W45 | 18 | 12456 | 8.26% | 34.91% |
1WW comparisons mean differential analysis between two isogenic counterparts in this research. For example, W12 means comparison between WT1 and WT2, and so forth.
2WG comparisons mean differential analysis between genetically modified rice varieties and their corresponding counterparts. For example, WG1 means comparison between WT1 and GM1, and so forth.
| C2 | C3 | C4 | C5* | |
|---|---|---|---|---|
| WG1 | 1 | 1 | 1 | 1 |
| WG2 | 2 | 3 | 4 | 3 |
| WG3 | 2 | 3 | 2 | 2 |
| WG4 | 2 | 3 | 4 | 5 |
| WG5 | 2 | 3 | 4 | 5 |
| WG6 | 2 | 2 | 3 | 4 |
| WG7 | 2 | 3 | 4 | 5 |
| WG8 | 2 | 3 | 4 | 5 |
*C2–C5 means clusters are classified into 2, 3, 4, and 5, respectively.
| Gene IDs | Times | Function | Compound IDs | Times | Metabolites |
|---|---|---|---|---|---|
| OS05G0144900 | 4 | Sucrose-phosphatase | C00022 | 3 | Pyruvate |
| OS06G0222100 | 4 | Trehalose-phosphate phosphatase 8 | C00025 | 3 | Glutamate |
| OS07G0463600 | 4 | Pentadecatrienyl resorcinol O-methyltransferase | C00072 | 3 | Ascorbate (vitamin C) |
| OS08G0459600 | 4 | 12-Oxophytodienoate reductase 7 | C00073 | 3 | Methionine |
| OS11G0135000 | 4 | Zinc induced facilitator | C00847 | 3 | Pyridoxate |
| OS10G0188100 | 4 | Transcription initiation factor IIF subunit alpha | C00864 | 3 | Pantothenate |
| OS02G0620400 | 4 | Uncharacterized | C01035 | 3 | 4-Guanidinobutanoate |
| OS02G0690000 | 4 | Uncharacterized | C01235 | 3 | Galactinol |
| C10646 | 3 | Lariciresinol |
| Description | Times | WG comparisons |
|---|---|---|
| GO overlap | ||
| Regulation of transcription, DNA dependent | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of primary metabolic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of nucleobase-containing compound metabolic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of nitrogen compound metabolic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of macromolecule metabolic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of cellular metabolic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| Regulation of biosynthetic process | 5 | WG2/WG4/WG5/WG7/WG8 |
| RNA biosynthetic process | 3 | WG4/WG6/WG8 |
| Defense response | 3 | WG1/WG2/WG3 |
| Regulation of biological process | 2 | WG2/WG8 |
| Programmed cell death | 2 | WG2/WG3 |
| Photosynthesis | 2 | WG1/WG7 |
| Lipid transport | 2 | WG2/WG7 |
| Death | 2 | WG2/WG3 |
| Cell death | 2 | WG2/WG3 |
| Apoptotic process | 2 | WG2/WG3 |
| KEGG overlap | ||
| Metabolic pathways | 8 | WG1–WG8 |
| Biosynthesis of secondary metabolites | 6 | WG1–WG6 |
| Starch and sucrose metabolism | 3 | WG2, WG5, WG8 |
| Ubiquitin mediated proteolysis | 2 | WG2, WG7 |
| Plant hormone signal transduction | 2 | WG2, WG7 |
| Photosynthesis | 2 | WG1, WG3 |
| WG comparisons | Gene ID | Description | Category |
|---|---|---|---|
| WG1 | OS10G0493600 OS07G0461900 | Alpha-galactosidase Acetylornithine aminotransferase | Metabolism/Defense response Metabolism |
| WG2 | OS06G0256500 | Glucose-6-phosphate | Metabolism |
| WG3 | OS12G0221700 OS07G0689600 | Leucine-rich repeat receptor Nicotianamine synthase | Signal reception and transduction Defense response |
| WG4 | OS03G0718000 | Anthranilate synthase beta subunit | Metabolism |
| WG5 | OS03G0136900 | Aconitate hydratase | Metabolism |
| WG6 | OS09G0567900 OS01G0316100 | Uridine nucleosidase Sarcosine oxidase | Metabolism Stress/Defense response |
| WG7 | OS03G0136900 | Aconitate hydratase | Metabolism |
| WG8 | – | – | – |
| GM materials | Negative control | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Stepsize | 3 | 5 | 10 | 20 | 30 | 50 | 3 | 5 | 10 | 20 | 30 | 50 |
| WG1 | 6* | 7 | 6 | 6 | 7 | 7 | 4 | 7 | 7 | 8 | 8 | 9 |
| WG2 | 3 | 7 | 6 | 7 | 4 | 6 | 7 | 9 | 9 | 6 | 6 | 6 |
| WG3 | 7 | 4 | 5 | 3 | 6 | 6 | 11 | 11 | 10 | 11 | 11 | 9 |
| WG4 | 3 | 5 | 2 | 5 | 6 | 6 | 11 | 12 | 10 | 10 | 11 | 9 |
| WG5 | 5 | 6 | 10 | 10 | 11 | 7 | 10 | 9 | 5 | 7 | 7 | 8 |
| WG6 | 7 | 6 | 6 | 9 | 9 | 8 | 9 | 12 | 9 | 11 | 10 | 7 |
| WG7 | 2 | 7 | 7 | 11 | 11 | 7 | 8 | 9 | 11 | 7 | 6 | 6 |
| WG8 | 12 | 10 | 6 | 10 | 9 | 8 | 5 | 7 | 8 | 7 | 6 | 7 |
*Number in bold indicates that this stepsize satisfied the condition that P value of the first cluster is the biggest in that of all the five clusters.
Supplemental Material
Download MS Word (566 KB)Related Research Data
Acknowledgments
We thank Prof. Wentao Xu, Prof. Jianbing Yan, and Prof. Daoxin Xie for providing useful suggestions.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Ethics Approval and Consent to Participate
Not applicable.
Availability of Data and Materials
RNA-seq read data have been deposited on the NCBI database under the accession number SRA230931-SRA230941.
Competing Interests
The authors declare that they have no competing interests.
Author Contributions
FW and WCG performed the experiment and analyzed the gene and metabolites mapping data. ZPY and XWJ performed the GO and KEGG analysis. ZYQ and ZSF conceived and designed the experiment. FW and ZSF provided the project resources. FW, WCG, and XWJ wrote the manuscript.
Supplementary material
Supplemental data for this article can be accessed on the publisher’s website.
Abbreviations
| ANOVA | analysis of variance |
| DEGs | differentially expressed genes |
| FPKM | fragments per kilobase million |
| GM | genetically modified |
| GO | gene ontology |
| KEGG | Kyoto Encyclopedia of Genes and Genomes |
| OPR | Oxophytodienoate reductase |
| PCA | principle component analysis |
| qPCR | quantitative polymerase chain reaction |
| TFIIF | transcription initiation factor IIF |
| WT | wild type |
| ZIF | zinc-induced facilitator |