Analysis of the Unintended Effects of the Bacillus thuringiensis Insecticidal Protein in Genetically Modiﬁed Rice Using Untargeted Transcriptomics

: The safety and unintended effects of genetically modiﬁed (GM) crops have been the focus of public attention. Transcriptome analysis is a powerful tool to assess the potential impact of genetic modiﬁcation on plant genomes. In this study, three transgenic (KMD, KF6, and TT51-1) and three non-transgenic (XS11, MH86, and MH63) rice varieties were assessed at the genomic and protein levels. The results of polymerase chain reaction (PCR) and Cry1Ab/1Ac speed test strips showed that the Bt gene was successfully expressed in transgenic rice. The results of RNA-seq analysis to analyze the unintended effects of transgenic Bt rice showed fewer differentially expressed genes (DEGs) between the transgenic and non-transgenic rice varieties than among the different varieties. Meanwhile, the results of principal component analysis and cluster analysis found no signiﬁcant genetic variation between the transgenic and non-transgenic rice varieties, except for the presence of Bt in transgenic rice. There were only two co-upregulated DEGs and no co-downregulated DEGs among three comparison groups. Although there were various DEGs among the groups, the two co-upregulated DEGs were not related to any signiﬁcantly enriched gene ontology (GO) term or Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, indicating that the differences among the subgroups were more likely caused by complex environmental or genetic factors, rather than unintended effects due to Bt expression. This study provides useful information to further explore the unexpected effects and safety of GM rice.


Introduction
Rice is the main food crop of almost half of the global population and has been historically important to the survival and development of human civilization [1].However, insect pests often cause huge losses to rice production and negatively impact food security [2,3], especially lepidopterans, such as the Asiatic rice borer (Chilo suppressalis), the yellow stem borer (Tryporyza incertulas), and the brown planthopper (Nilaparvata lugens).
Chemical pesticides are conventionally applied to control pests of rice crops, but these products are environmental pollutants that have been linked to drug resistance, destruction of ecosystems, and rampant pest re-infestation.Therefore, various strategies have been proposed to improve rice production without the use of chemical pesticides [4].
Transgenic technologies have been applied as cost-effective methods to increase crop production by improving resistance to insect pests and diseases, tolerance to herbicides, and nutritional value [5,6].An insecticidal crystal protein derived from the Bt protoxin produced by Bacillus thuringiensis exhibits highly specific insecticidal activity and is commonly used in research of transgenic insect-resistant rice [7].The crystal protein derived from the Bt protoxin penetrates the intestine of insect pests.From 1996 to 2019, Bt maize was the first large-scale genetically modified (GM) crop and is currently grown on an estimated 190.4 million hectares in more than 70 countries worldwide [8,9].However, the rapid increase in cultivation of GM crops has led to increasing global concern about food safety [10].It is generally agreed that biosafety assessment of GM crops and related products should also include potential long-term unintended effects, such as genomic instability and random integration of transgenes [11][12][13].Therefore, before planting GM crops, such as Bt rice, risk evaluations should be conducted to assess potential harmful effects on agroecosystems [14,15].
Genomic characterization is typically employed to assess the safety of GM crops.The development of omics techniques, including transcriptomics [6,[16][17][18], proteomics [19,20], and metabolomics [21][22][23], has facilitated comprehensive analyses of the expected and unintended effects of transgenic events in GM crops at the transcript, protein, and metabolite levels [24].Based on this information, the expected and unintended effects of transgenes can be easily identified [25,26], and risk assessments can be conducted in a timely manner to assess the safety of GM crops and related products.
Transcriptomics is the systematic study of transcriptional profiles, which aims to elucidate the molecular mechanisms of complex biological pathways and trait regulatory networks.Transcriptomics not only detects transcripts corresponding to existing genome sequences, but also discovers and quantifies new transcripts, making it more advantageous for the study of selective splicing events, new genes, and transcripts.Currently, commonly used techniques for transcriptomics sequencing analysis include hybridizationbased techniques, such as DNA microarray [27]; label-based techniques, such as serial analysis of gene expression (SAGE) [28]; and direct sequencing-based techniques, such as RNA sequencing (RNA-seq) [29].DNA microarray technology is mature and has a large database, but it is only suitable for detecting known sequences and has low sensitivity; SAGE technology has a high demand for mRNA, while RNA-seq does not require known sequences and is able to detect the overall transcriptional activity of any species at the single nucleotide level [30].In recent years, RNA sequencing (RNA-seq) analysis has been widely applied in omics studies of the unintended expression of transgenes.For example, Peng et al. used RNA-seq analysis to compare the miRNA profiles of transgenic Bt and 5-enolpyruvylshikimate 3-phosphate synthase (EPSPS) varieties and non-transgenic rice varieties [31].Using transcriptome analysis, Guo et al. revealed differences in the response and tolerance mechanisms of EPSPS and GAT in transgenic soybeans [32].Moreover, a study by Wang et al. using RNA-seq analysis found fewer DEGs between transgenic and non-transgenic maize than among the different varieties [24].
The aim of the present study was to detect Bt protein levels in samples with the use of Bt-Cry1Ab/1Ac detection test strips.Specific primers were designed for three transgenic rice varieties (KMD, KF6, and TT51-1) and three non-transgenic rice varieties (XS11, MH86, and MH63).Total RNA of transgenic Bt rice was collected for RNA-seq analysis to identify DEGs.The biological processes, cellular locations, and molecular functions of the DEGs were determined through gene ontology (GO) analysis in reference to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.An experimental flow chart is shown in Figure 1.The results of this study provide valuable references for genomics-based evaluation of the safety of GM rice.
The rice seeds were selected with full shape and soaked in 70% alcohol for 1 min for surface disinfection.The seeds were then soaked in 2% sodium hypochlorite solution for 15 min and rinsed with water 4 times.The seeds were placed in an incubator at 28 °C for germination.The better germinated seeds were selected and sown in 3 plastic pots (10 seeds per pot) for each variety (three transgenic varieties and three non-transgenic rice varieties) and placed in the incubator.The incubation conditions were set to 28 ℃/16 h light, 24 ℃/8 h dark, and 80% humidity.Thirty days later, each sample was pooled from the leaves collected from three plastic pots (each pot needs to be collected) for each variety.The same sample was taken for DNA, RNA, and protein extraction, respectively.

DNA Extraction, PCR Analysis, and Test Strips
DNA was extracted from rice samples using Plant Genomic DNA Extraction Kit (Tiangen Biochemical Technology Co., Ltd.Beijing, China), and the purity of DNA was determined using the ratio of the absorbance at 260 and 280 nm using a spectrophotometer (Nanodrop 2000, ThermoFisher Scientific, Waltham, MA, USA).The DNA samples were amplified using qualitative polymerase chain reaction (PCR).The 25 uL PCR amplification system comprised 10 × PCR buffer (Mg 2+ Plus) 2.5 µL, dNTP mixture 2 µL, 0.5 µL each of forward and reverse primers (10 µmol/L) at the final concentration of 200 nmol/L, rTaq DNA polymerase (5 U/µL) 0.15 µL, and DNA solution 2 µL; the reaction mixture was supplemented with ddH2O to make up the volume to 25 µL.The amplification program for each primer set is listed in Supplementary Table S1.
Target proteins were detected by test strips.Leaves were placed in a 1.5 mL centrifuge tube, the leaves were triturated, 1 mL buffer was added, and a test strip (Shanghai
The rice seeds were selected with full shape and soaked in 70% alcohol for 1 min for surface disinfection.The seeds were then soaked in 2% sodium hypochlorite solution for 15 min and rinsed with water 4 times.The seeds were placed in an incubator at 28 • C for germination.The better germinated seeds were selected and sown in 3 plastic pots (10 seeds per pot) for each variety (three transgenic varieties and three non-transgenic rice varieties) and placed in the incubator.The incubation conditions were set to 28 °C/16 h light, 24 °C/8 h dark, and 80% humidity.Thirty days later, each sample was pooled from the leaves collected from three plastic pots (each pot needs to be collected) for each variety.The same sample was taken for DNA, RNA, and protein extraction, respectively.

DNA Extraction, PCR Analysis, and Test Strips
DNA was extracted from rice samples using Plant Genomic DNA Extraction Kit (Tiangen Biochemical Technology Co., Ltd.Beijing, China), and the purity of DNA was determined using the ratio of the absorbance at 260 and 280 nm using a spectrophotometer (Nanodrop 2000, ThermoFisher Scientific, Waltham, MA, USA).The DNA samples were amplified using qualitative polymerase chain reaction (PCR).The 25 uL PCR amplification system comprised 10 × PCR buffer (Mg 2+ Plus) 2.5 µL, dNTP mixture 2 µL, 0.5 µL each of forward and reverse primers (10 µmol/L) at the final concentration of 200 nmol/L, rTaq DNA polymerase (5 U/µL) 0.15 µL, and DNA solution 2 µL; the reaction mixture was supplemented with ddH 2 O to make up the volume to 25 µL.The amplification program for each primer set is listed in Supplementary Table S1.
Target proteins were detected by test strips.Leaves were placed in a 1.5 mL centrifuge tube, the leaves were triturated, 1 mL buffer was added, and a test strip (Shanghai Youlong Biotechnology Co., Ltd.Shanghai, China) was inserted into the sample mixture; observe the results after 5 min.

RNA Extraction, Library Preparation, and Sequencing
Total RNA from rice leaves was extracted using an RNAprep pure Plant kit (TIANGEN Biotech, Beijing, China), and digested at 37 • C for 1 h using DNaseI (Ambion, Austin, TX, USA); the purified RNA was assayed for RNA quality using Qubit ® 2.0.The quality of the extracted RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).Library preparation with Illumina sequencing was performed by following the protocol plan.Extracted total RNA was subjected to mRNA enrichment, followed by short fragmentation of the enriched mRNA, cDNA first-strand synthesis with random primers, followed by cDNA second-strand synthesis, end repair, and addition of poly A after purification using a QiaQuick PCR kit (Qiagen, Duesseldorf, Germany), and ligation to Illumina sequencing adapters.The ligated products were subjected to agarose gel electrophoresis, PCR amplification, and sequencing using Illumina HiSeqTM4000.

RNA-Seq Data Analysis
The rice genome was aligned with the reference genome using bowtie2.Differential expression analysis of XS11 vs. KMD, MH86 vs. KF6, and MH63 vs. TT51-1 data was performed using the software GFOLD [38], the screening conditions for differential expression genes were set at |GFOLD_value| > 1.
GO is an international standardized gene function classification system, which is a database applicable to each species to qualify and describe the functions of genes and proteins.Using the GO database, genes are classified according to the biological processes, cellular components, and molecular functions.GO enrichment analysis of DEGs was performed using the R package (clusterProfiler, Vienna, Austria), and GO terms with p < 0.05 were considered significantly enriched in DEGs [39,40].
KEGG is an integrated database resource for biological interpretation of genomic sequences and other high-throughput data.Molecular functions of genes and proteins are associated with immediate homology groups and stored in the KEGG Orthology (KO) database.The R package (clusterProfiler) was used to detect the statistical enrichment of DEGs in the KEGG pathway.KEGG pathways with p < 0.05 were considered significantly enriched in DEGs [40,41].

Bt Protein Expression in GM Rice
The synthetic Cry1Ab/Ac gene encodes an insecticidal protein of B. thuringiensis and is the most widely used exogenous protein for insecticide genetic engineering worldwide.The Bt crystal protein is a δ-endotoxin derived from the Bt protoxin produced by B. thuringiensis in the early stage of endospore formation.Diagrams of the vectors encoding the transgenes of the three transgenic rice varieties are shown in Figure 2A-C.Primers to detect exogenous genes were based on the Cry1Ab/Ac sequence.The rice sucrose phosphate synthase (SPS) gene was used as an internal reference.The SPS sequences of all transgenic and nontransgenic (control) samples were amplified by PCR.Target fragments of the exogenous genes of the three transgenic rice varieties were successfully amplified, but not those of the non-transgenic rice.Further, all six samples were examined using Cry1Ab/Ac test strips.Notably, the test strips of transgenic rice exhibited two red lines, while those of the non-transgenic rice had only one red line, indicating Bt protein expression by the three transgenic varieties (Figure 2D).The relative RNA expression of Bt obtained from RNA-seq is shown in Supplementary Figure S1.

Evaluation of Rice Transcriptome Sequencing Data
Transcriptome sequencing of rice leaves was performed to identify differen gene expression patterns between transgenic (KMD, KF6, and TT51-1) non-transgenic (XS11, MH86, and MH63) rice varieties.RNA-seq analysis of fi quality data of the six samples yielded 12112186-12390376 raw reads with a GC co of 47-52%.As shown in Table 1, the clean reads of each sample shared similar 83.24-87.06%with the Oryza sativa Japonica Group (Japanese rice) genome asse IRGSP-1.0 from National Institute of Agrobiological Sciences (GCA_001433 GCF_001433935.1).The gene expression levels of the rice samples are expressed as ments per kilobase of transcript per million mapped fragments (FPKM).As show Supplementary Figure S2, gene expression levels were similar among the six sam The measured data were accurate with good quality and sufficient for further analy

Evaluation of Rice Transcriptome Sequencing Data
Transcriptome sequencing of rice leaves was performed to identify difference in gene expression patterns between transgenic (KMD, KF6, and TT51-1) and non-transgenic (XS11, MH86, and MH63) rice varieties.RNA-seq analysis of filtered quality data of the six samples yielded 12112186-12390376 raw reads with a GC content of 47-52%.As shown in Table 1, the clean reads of each sample shared similarity of 83.24-87.06%with the Oryza sativa Japonica Group (Japanese rice) genome assembly IRGSP-1.0 from National Institute of Agrobiological Sciences (GCA_001433935.1 GCF_001433935.1).The gene expression levels of the rice samples are expressed as fragments per kilobase of transcript per million mapped fragments (FPKM).As shown in Supplementary Figure S2, gene expression levels were similar among the six samples.The measured data were accurate with good quality and sufficient for further analysis.

DEGs between the Transgenic and Non-Transgenic Rice Varieties
Principal component analysis (PCA) was performed to assess the genetic variation between the transgenic and non-transgenic rice varieties.Three-dimensional plots of the PCA results of the different rice varieties are presented in Figure 3. Notably, XS11 and KMD were significantly separated from MH86, KF6, MH63, and TT51-1 (PC1 axis), accounting for 72.70% of the variability among the samples.Notably, XS11 and KMD are Japonica varieties, while MH86, KF6, MH63, and TT51-1 are Indica varieties.As shown in Figure 3, MH86 clustered with KF6 and MH63 with TT51-1 along the PC2 axis.In addition, XS11 and KMD, which accounted for 12.04% of the variability among the samples, were clustered together at the most negative range of the PC2 axis.Through combined analysis, there was no significant difference between the transgenic and non-transgenic (control) groups.The clustering plot in Supplementary Figure S3A shows clear differences among the samples.Although all six samples belonged to the same category, XS11 was closely related to KMD, while MH86, KF6, MH63, and TT51-1 were closely related to one another.As shown in Supplementary Figure S3B, the Pearson's correlation coefficient(r) was 0.94 for XS11 and KMD, 0.98 for MH86 and KF6, and 0.95 for MH86 and TT51-1, indicating no significant difference between the transgenic and non-transgenic (control) groups.
Processes 2023, 9, x FOR PEER REVIEW 6 o

DEGs between the Transgenic and Non-Transgenic Rice Varieties
Principal component analysis (PCA) was performed to assess the genetic variat between the transgenic and non-transgenic rice varieties.Three-dimensional plots of PCA results of the different rice varieties are presented in Figure 3. Notably, XS11 a KMD were significantly separated from MH86, KF6, MH63, and TT51-1 (PC1 axis), counting for 72.70% of the variability among the samples.Notably, XS11 and KMD Japonica varieties, while MH86, KF6, MH63, and TT51-1 are Indica varieties.As shown Figure 3, MH86 clustered with KF6 and MH63 with TT51-1 along the PC2 axis.In ad tion, XS11 and KMD, which accounted for 12.04% of the variability among the samp were clustered together at the most negative range of the PC2 axis.Through combin analysis, there was no significant difference between the transgenic and non-transge (control) groups.The clustering plot in Supplementary Figure S3A shows clear dif ences among the samples.Although all six samples belonged to the same category, X was closely related to KMD, while MH86, KF6, MH63, and TT51-1 were closely related one another.As shown in Supplementary Figure S3B, the Pearson's correlation coe cient(r) was 0.94 for XS11 and KMD, 0.98 for MH86 and KF6, and 0.95 for MH86 a TT51-1, indicating no significant difference between the transgenic and non-transge (control) groups.
(Figure 4D).Os04g0405100 was annotated as a SAM domain family protein and Os06g0115300 as acyl coenzyme A binding.In addition, 3351 DEGs were identified between XS11 and MH86, while 2893 DEGs were identified between XS11 and MH63.The number of DEGs was lower between transgenic and non-transgenic rice varieties than among the different varieties.

GO and KEGG Analysis of DEGs
GO describes three distinct functional classes of gene products: molecular function, cellular component, and biological process.GO and KEGG enrichment analysis was performed to explore the biological functions of the DEGs of transgenic and non-transgenic rice varieties.The two co-upregulated DEGs were not associated with any GO terms.The biological functions of DEGs for XS11 vs. KMD included "cellular macromolecular biosynthetic processes" and "small molecule metabolic processes".The biological functions of DEGs for MH86 vs. KF6 included "diterpenoid metabolic processes" and "terpenoid metabolic processes".The biological functions of DEGs for MH63 vs. TT51-1 included "GTP binding" and "guanyl nucleotide binding" (Figure 5).Notably, there were no common significantly enriched GO terms among the three pairs of rice varieties.

GO and KEGG Analysis of DEGs
GO describes three distinct functional classes of gene products: molecular function, cellular component, and biological process.GO and KEGG enrichment analysis was performed to explore the biological functions of the DEGs of transgenic and non-transgenic rice varieties.The two co-upregulated DEGs were not associated with any GO terms.The biological functions of DEGs for XS11 vs. KMD included "cellular macromolecular biosynthetic processes" and "small molecule metabolic processes".The biological functions of DEGs for MH86 vs. KF6 included "diterpenoid metabolic processes" and "terpenoid metabolic processes".The biological functions of DEGs for MH63 vs. TT51-1 included "GTP binding" and "guanyl nucleotide binding" (Figure 5).Notably, there were no common significantly enriched GO terms among the three pairs of rice varieties.KEGG pathways associated with the DEGs are shown in Figure 6.The two co-upregulated genes were not significantly enriched in any KEGG pathway.The DEGs for XS11 vs. KMD were significantly enriched in 13 KEGG pathways, which included "genetic information processing", "metabolism", "ribosome", "biosynthesis of cofactor", and "photosynthesis".The DEGs for MH86 vs. KF6 were significantly enriched in 10 KEGG pathways, which included "cellular processes", "genetic information processing", "metabolism", "diterpenoid biosynthesis", "motor protein", and "DNA replication".The DEGs for MH63 vs. TT51-1 were significantly enriched in 10 KEGG pathways, which included "metabolism", "genetic information processing", "biotin metabolism", "glutathione metabolism", and "linoleic acid metabolism".These pathways are involved in basic plant metabolism.None of the significantly enriched pathways were shared among the three pairs of rice varieties.KEGG pathways associated with the DEGs are shown in Figure 6.The two coupregulated genes were not significantly enriched in any KEGG pathway.The DEGs for XS11 vs. KMD were significantly enriched in 13 KEGG pathways, which included "genetic information processing", "metabolism", "ribosome", "biosynthesis of cofactor", and "photosynthesis".The DEGs for MH86 vs. KF6 were significantly enriched in 10 KEGG pathways, which included "cellular processes", "genetic information processing", "metabolism", "diterpenoid biosynthesis", "motor protein", and "DNA replication".The DEGs for MH63 vs. TT51-1 were significantly enriched in 10 KEGG pathways, which included "metabolism", "genetic information processing", "biotin metabolism", "glutathione metabolism", and "linoleic acid metabolism".These pathways are involved in basic plant metabolism.None of the significantly enriched pathways were shared among the three pairs of rice varieties.

Discussion
Transgenic technology is increasingly used to improve crop quality and yield.At present, GM crops are farmed in more than 70 countries worldwide.However, the safety of GM crops and related products remains a major concern.Hence, the risks of GM crops to human health and potential environmental impacts must be carefully assessed.Nonetheless, it is difficult to predict the unintended effects of transgenes [24][25][26], such as genetic disruptions leading to sequence changes and the production of new proteins or metabolites, which could impair the safety of GM crops [12,42].To address these issues, omics technologies, such as transcriptomics, proteomics, and metabolomics, can facilitate the analysis of the unintended effects of GM crops [43][44][45][46].The emergence of Omics provides a new method for evaluating the safety of GM foods, and transcriptomics can be used for low-abundance transcripts, deep exploration of new genes, identification of gene families, determination of metabolic pathways, and other analyses, which have been widely used in many fields [47,48].
Omics technologies generate large volumes of data, and thus, PCA was conducted to simplify presentation of relevant information.Here, three transgenic (KMD, KF6, and TT51-1) and three non-transgenic (XS11, MH86, and MH63) rice varieties were analyzed to identify potential unintended effects of GM rice.The selected samples included Japonica (XS11 and KMD) and Indica (MH68, KF6, MH63, and TT51-1) rice varieties.The PCA results showed similarity of the transgenic and non-transgenic rice varieties, while the results of cluster analysis demonstrated the separation of the Japonica and Indica varieties, as XS11 clustered with KMD, while MH86, KF6, MH63, and TT51-1 clustered together, consistent with expectations and in agreement with the results of a previous study [3].Notably, there was no significant separation of the transgenic and non-transgenic Indica or Japonica rice varieties, in accordance with previous studies of GM barley [49].
Overall, there were fewer DEGs between the transgenic and non-transgenic (XS11 vs. KMD, MH86 vs. KF6, and MH63 vs. TT51-1) than between the Japonica and Indica rice varieties (XS11 vs. MH86 and XS11 vs. MH63), similarly to the results of previous studies of GM maize, rice and soybean [24,43,50].Another study found more DEGs between rice roots and leaves than between transgenic and non-transgenic rice varieties [51].The Os04g0405100 and Os06g0115300 (OsACBP2) genes, which were co-upregulated in XS11 vs. KMD, MH86 vs. KF6, and MH63 vs. TT51-1, were annotated to SAM domain family protein and acyl coenzyme A binding, respectively.The evolutionary conserved and ubiquitous family of acyl-CoA-binding proteins (ACBPs) are involved in basal cellular functions [52].ACBP1 and ACBP2 of Arabidopsis thaliana play important roles in phospholipid membrane repair in response to heavy metal stress and are critical to early embryonic development [53][54][55][56].Unlike ACBPs of members of the Arabidopsis genus, which are involved in different pathways in response to both biotic and abiotic stressors, ACBP2 of Oryza sativa is not induced by stressors [57].
The two co-upregulated genes were not associated with enriched GO terms or KEGG pathways.Some of the DEGs identified in this study were enriched in various biological processes, including cellular macromolecular biosynthesis, small molecule metabolic processes, carbohydrate derivative metabolic processes, diterpene metabolic processes, GTP binding, and among others.The DEGs of XS11 and KMD were significantly enriched in various KEGG pathways, including ribosome and cofactor biosynthesis, while those of MH86 and KF6 were significantly enriched in motor protein and DNA replication pathways, and the DEGs of MH63 and TT51-1 were significantly enriched in biotin metabolism and glutathione metabolism, which are basic metabolic pathways in plants.These data on the profiling of transgenic rice reveal that some differences exist compared to non-transgenic rice.This may be due to the fact that the selection process for transgenic rice is based on phenotypic and compositional equivalence with a close comparator line, and then the introduction of the new trait into a superior line by a number of crosses, rather than the proper expression of the new trait [58].Some environmental factors and insertional effects of transgenic genes have also been shown to have a greater influence than transgenic gene expression [59,60].In the future, we plan to sow transgenic and non-transgenic crops in different environments at different years to reveal the environmental and genetic effects.
In summary, transcriptomics is an effective method to assess the unexpected effects of GM crops.In the present study, RNA-seq analysis was conducted to analyze gene expression levels in transgenic vs. non-transgenic rice varieties.Previous studies have reported that the expression levels of DEGs in transgenic vs. non-transgenic rice were within expected ranges as compared to inter-varietal differences.Moreover, although attempts were made to exclude environmental factors as much as possible, other factors, such as natural variation and genetic background, could have influenced the differences in insertional expression of Bt.Overall, trans-BT GM rice had no widespread unintended effects.The results of this study provide new perspectives for safety evaluation of GM rice.

Figure 3 .
Figure 3. Principal component analysis of samples.PCA results showing similarity of the tr genic and non-transgenic rice varieties.The Japonica varieties were separated from the Indica va ties.

Figure 3 .
Figure 3. Principal component analysis of samples.PCA results showing similarity of the transgenic and non-transgenic rice varieties.The Japonica varieties were separated from the Indica varieties.

Figure 5 .
Figure 5. GO analysis of DEGs.The significantly enriched GO terms of the DEGs among the three groups are shown (p < 0.05, from small to large).The size of the dots represents the number of genes.Biological processes are indicated by green lines and molecular function by orange lines.

Figure 5 .
Figure 5. GO analysis of DEGs.The significantly enriched GO terms of the DEGs among the three groups are shown (p < 0.05, from small to large).The size of the dots represents the number of genes.Biological processes are indicated by green lines and molecular function by orange lines.

Table 1 .
RNA-sequencing and transcriptome mapping results.

Table 1 .
RNA-sequencing and transcriptome mapping results.