Analysis of transcriptomic differences between NK603 maize and near-isogenic varieties using RNA sequencing and RT-qPCR

The insertion of a transgene into a plant organism can, in addition to the intended effects, lead to unintended effects in the plants. To uncover such effects, we compared maize grains of two genetically modified varieties containing NK603 (AG8025RR2, AG9045RR2) to their non-transgenic counterparts (AG8025conv, AG9045conv) using high-throughput RNA sequencing. Moreover, in-depth analysis of these data was performed to reveal the biological meaning of detected differences. Uniquely mapped reads corresponded to 29,146 and 33,420 counts in the AG8025 and AG9045 varieties, respectively. An analysis using the R-Bioconductor package EdgeR revealed 3534 and 694 DEGs (significant differentially expressed genes) between the varieties AG8025RR2 and AG9045RR2, respectively, and their non-transgenic counterparts. Furthermore, a Deseq2 package revealed 2477 and 440 DEGs between AG8025RR2 and AG9045RR2, respectively, and their counterparts. We were able to confirm the RNA-seq results by the analysis of two randomly selected genes using RT-qPCR (reverse transcription quantitative PCR). PCA and heatmap analysis confirmed a robust data set that differentiates the genotypes even by transgenic event. A detailed analysis of the DEGs was performed by the functional annotation of GO (Gene Ontology), annotation/enrichment analysis of KEGG (Kyoto Encyclopedia of Genes and Genomes) ontologies and functional classification of resulting key genes using the DAVID Bioinformatics Package. Several biological processes and metabolic pathways were found to be significantly different in both variety pairs. Overall, our data clearly demonstrate substantial differences between the analyzed transgenic varieties and their non-transgenic counterparts. These differences indicate that several unintended effects have occurred as a result of NK603 integration. Heatmap data imply that most of the transgenic insert effects are variety-dependent. However, identified key genes involved in affected pathways of both variety pairs show that transgenic independent effects cannot be excluded. Further research of different NK603 varieties is necessary to clarify the role of internal and external influences on gene expression. Nevertheless, our study suggests that RNA-seq analysis can be utilized as a tool to characterize unintended genetic effects in transgenic plants and may also be useful in the safety assessment and authorization of genetically modified (GM) plants.

Roundup Ready maize (NK603) was analyzed. The NK603 transgene consists of two cassettes, which were inserted by microparticle bombardment. Both cassettes contain a 5-enolpyruvylshikimate-3-phosphate synthase (EPSPS) gene obtained from the soil bacterium Agrobacterium sp. strain CP4. In the first cassette, the gene is regulated by the rice (Oryza sativa) actin promoter. The second cassette has nearly the same composition as the first one, but it is regulated by an enhanced 35S promoter obtained from the cauliflower mosaic virus. The CP4-EPSPS transgene confers plant tolerance to the herbicide glyphosate, commercially sold as Roundup ® . The function of the EPSPS protein, its toxicity and allergenicity as well as metabolites have been examined in several projects [9,29,39,70]. In the field of genetic engineering, there is still a lack of clarity about transgene integration mechanisms and about which biological processes take place in the plant after successful transgene integration [46,74]. Since DNA integration can lead to gene disruptions, DNA rearrangements, or the production of new proteins [71], the insertion of a transgene into a plant organism can, in addition to the intended effects (e.g., herbicide resistance), lead to unintended effects in the plants [38,50,74]. Such unintended effects need to be identified and evaluated to detect and prevent possible adverse effects [28].
For the authorization and approval of GM crops, it is crucial to avoid the occurrence of unintended effects. Thus, investigations of GM crops on a molecular biology level are essential to enable a deeper understanding of plant gene interactions. Currently, unintended effects are primarily investigated by targeted approaches, in which GM crops are compared to their near-isogenic non-GM counterparts. Molecular, compositional, phenotypic, and agronomic analyses are performed in order to identify similarities and differences between the crops. Significant differences between a GM crop and its conventional counterpart indicate an unintended effect, which requires further investigation [30]. Targeted approaches have the disadvantage that important differences can be overlooked, as only selected features or chemical/nutritional compounds are investigated. Untargeted omics approaches, e.g., the analysis of gene expression in its entirety, enable the identification of broad coordinated trends that cannot be discerned by targeted approaches [54]. Thus, new non-targeted profiling approaches are an appropriate tool to complement current investigations regarding unintended effects of transgenic plants [31,63,71].
There are several studies investigating unintended effects in GM crops using different omics approaches and microarrays. In these studies, it was found that transgenic inserts of different species of plants might affect the overall expression of other endogenous genes [1,2,7,8,12,17,19,32,35,40,45,49,60,84,90]. However, several authors of these studies mention that the changes in gene expression and protein distribution caused by genetic modification were smaller than those caused by environmental factors or natural variations. In addition, no differences were found in some of these studies. For instance, [21] and 2009, compared grown MON810 maize varieties with some comparable varieties using microarrays under in vitro conditions and conventional agricultural field conditions. In these studies, no gene was found to be significantly differentially expressed in any of the variety pairs tested under in vitro or conventional field conditions [20,21].
Differences in lignin content between MON810 maize and comparators were reported by Saxena and Stotzky as well as Poerschmann [65,73]. Herrero et al. detected significant differences in enantiomeric amino acid composition between two MON810 varieties and comparators but not in a third tested variety pair (MON810 vs. non-GM) [41]. Agapito et al. conducted a comparative analysis of MON810 maize and comparators grown under different agroecosystem conditions in Brazil using two-dimensional gel electrophoresis combined with mass spectrometry, which identified 32 differentially expressed proteins [1]. La Paz et al. performed RNA sequencing on MON810 maize and comparators, which revealed 140 differentially expressed genes. A slight, but significant, delay in seed and plant maturation of MON810 maize plants was also observed and, thus, used as an explanation for the detected differences [49]. [36].
Unintended effects in NK603 maize were described by several authors [2,8,12,37,40,58]. Agapito-Tenfen et al. analyzed two maize events, NK603 and MON89034, as single events and as stacked varieties. Proteomics was used for the identification of unintended effects. Twentytwo proteins were detected to be differentially expressed in stacked and single GM events compared to a nearisogenic non-GM maize and to a landrace variety [2]. Mesnage et al. generated proteome profiles of Roundup Ready maize and near-isogenic maize varieties using molecular profiling. A comparison of GM maize vs. the near-isogenic maize revealed alterations in the levels of enzymes of the glycolysis and TCA cycle pathways, which can be interpreted as an imbalance in energy metabolism [58]. Benevenuto et al. performed proteomic and metabolomic analyses on NK603 maize plants under abiotic stress conditions. Twenty differentially modulated proteins were identified between GM and non-GM hybrids under water deficiency conditions and herbicide sprays [12]. Barros et al. compared two Roundup Ready maize varieties with a near-isogenic non-GM variety using transcriptomics, proteomics, and metabolomics. The results showed that the environment had a major effect on protein and gene expression and metabolite production [8]. Harrigan et al. analyzed differences in the metabolome between NK603 hybrids, corresponding negative segregants, and conventional comparator hybrids. It was concluded that the largest effects on metabolomic variation were due to growing location and genomic differences associated with backcrossing practices [37].
In previous studies, we identified a silent mutation in the coding region of the cry1ab gene in different stacked MON810 maize varieties using high resolution melting (HRM) analysis, Sanger sequencing, and amplicon sequencing [10,11]. We have also performed a molecular analysis of NK603 maize and identified two insertions, which are not present in the NK603 patent sequence. These insertions were located in the transgenic promoter region; therefore, they may have an effect on the promoter activity and, consequently, also on transgene expression [14]. However, our analysis of the 5´-end of MON810 as a single event with Scorpion probes revealed no unintended effects or mutations [62].
In this study, we focused on a non-targeted transcriptomics approach (RNA-sequencing) for the identification of possible unintended effects in herbicide-resistant maize (NK603 maize). To the best of our knowledge, this is the first study in which NK603 maize varieties and corresponding non-GM counterparts were analysed by RNA sequencing. We compared the entire transcriptome of RR maize with near-isogenic varieties to characterize gene candidates that may be differentially expressed between RR maize and near-isogenic varieties in two different genetic backgrounds. We conducted a principal component analysis (PCA) as well as a heatmap to evaluate the variance structure of our data. Moreover, we performed GO, KEGG annotations/enrichment and DAVID analyses, which can identify biological functions and metabolic pathways particularly affected by the modification of gene expression. Our results indicate several unintended effects that may have occurred as a result of transgene integration.

Rationale and experimental design
The main objective of this study was to compare gene expression profiles of two different NK603 varieties (AG8025RR and AG9045RR) with the profiles of corresponding near-isogenic varieties that do not contain any transgenes (AG8025conv and AG9045conv, respectively) and investigate whether there are consistent alterations in the NK603 varieties. For this purpose, we used high-throughput RNA sequencing as well as RT-qPCR approach. The pipeline used for the analysis of this study is depicted in Fig. 1. DGE analysis was performed by two different statistical packages (DeSeq2, EdgeR).
Further, we performed a singular enrichment analysis (SEA) to identify if special gene classes and gene interacting networks that are overrepresented among the differentially expressed genes, and we identified the biological functions of differentially expressed genes using GO annotations. In order to search for shared GO terms in the two NK603 genetic backgrounds, we have performed a second joint analysis using REViGO tool followed by a functional pathway analysis with KEGG annotations. A KEGG metabolic pathways analysis was also conducted. These analyses were performed for each variety pair (NK603 and conventional counterpart) separately. Finally, we compared the genes of the overrepresented gene classes (SEA) with the genes involved in significantly affected metabolic pathways (KEGG) and conducted a DAVID analysis with these key genes to clarify the relationship between these genes.
Next to this, we performed RT-qPCR for determining if some of the detected differences in gene expression of AG8025RR2 vs. AG8025conv were also present in AG9045RR2 vs. AG9045conv. In parallel to these analyses, we have performed a PCA as well as a heatmap to evaluate the variance structure of our data set.  insert which enables glyphosate herbicide tolerance) and the near-isogenic varieties AG8025 (conventional counterpart of AG8025RR2) and AG9045 (conventional counterpart of AG9045RR2) were obtained from the Brazilian market (Sementes Agroceres). AG8025RR2 and its conventional counterpart AG8025 have the same genetic background since they are produced from the same endogamic parental lines. The same is for AG9045RR2 and its conventional counterpart AG9045. The AG8025 variety is the hybrid progeny of the single-cross between maternal endogamous line "A" with the paternal endogamous line "B". Thus, the tested hybrid variety grains have high genetic similarity (AB genotype). The AG9045 variety is the hybrid progeny of the single-cross between maternal endogamous line "C" with the paternal endogamous line "D" resulting in CD genotype. Presence/absence of the transgenic NK603 insert was confirmed by qPCR in both, GM and conventional maize grains. The obtained maize hybrids were directly used for the transcriptome analysis.

Plant material
We have not received any information about the exact place of cultivation, type of soil, the possible application of fertilizers or pesticides. However, it is confirmed that these varieties were produced in Brazil in 2012.

RNA extraction
Total RNA from pools of ten maize grains (of each variety) was isolated based on protocols of Barros et al. and Cheng et al. [8,15]. Maize grains were homogenized in liquid nitrogen using mortar and pestle. 500 mg of the homogenized grains was transferred into a 15 ml tube and 5 ml of 60 °C warm extraction buffer (2% CTAB, 2% PVP), 2 M NaCl, 0.9 mM DEPC, 0.5 mM spermidine, 100 mM Tris [pH = 8]) and 100 µl of mercaptoethanol were added. After vortexing, 5 ml of chloroform-isoamylalcohol (24:1; CIA) was added. The mixture was vortexed and centrifuged (

RNA sequencing and data analysis
Three pools of variety AG8025conv and three pools of variety AG8025RR2 were subjected to standard Illumina library preparation using the NEBNext Ultra RNA library prep kit according to the manufacturer's instructions. Three pools of varieties 9045conv and 9045RR2, respectively were treated in the same way. The resulting six cDNA libraries of each conv, RR2 variety pair were paired end sequenced (125 bp) using an Illumina HiSeq2500 machine at the Vienna Biocenter Core Facilities NGS unit (https ://www.vbcf.ac.at). Reads, which passed basic quality control (Illumina chastity filter), were preprocessed: Adapters were removed with cutadapt (https ://cutad apt. readt hedoc s.io/en/stabl e/guide .html) and the reads were filtered against the rDNA using Bowtie2 (https ://bowti e-bio.sourc eforg e.net/bowti e2/index .shtml ), which is a very sensitive contaminants database. The remaining reads were paired end aligned with STAR [25] against the B73 reference genome of Zea mays (AGPv3).

Statistical analyses
The Principal Component Analysis was performed using 'stats and 'ggplot2′ libraries while the heatmap clustering analysis was conducted using 'pheatmap' library in R environment aiming to find gene expression patterns across the different varieties. All clean read counts for each sample (n = 12) were used in these analyses. In order to determine DEGs, the two NK603 varieties AG9045RR2 and AG8025RR2 were compared to their isogenic counterparts AG9045conv and AG8025conv, respectively. Data were calculated by two tests, DEseq2 and EdgeR, using the software packages Bioconductor and Galaxy. These tests are among the best and most used performance tools for RNA-seq analysis with low numbers of false positives and reliable gene-wise dispersion estimates across all samples [4,53,72]. DEseq2 and EdgeR analysis are based on the assumption that the data follow a negative binomial distribution [27]. Using the raw counts, the data were normalized and transformed to correct for dispersion artifacts and variability within the compared groups or to account for differences in sequencing depths. Genes without any counts were removed. As significance level an adjusted p-value controlled for multiple testing using Benjamini and Hochberg's correction-with a false discovery rate (FDR) below or equal to ≤ 0.05 was taken for the characterization of DEGs [13,67,75,76]. In addition, a general cut-off threshold of a log2-fold change (log2 FC) ≥ + 1 as well as a log2 FC ≤ -1 was used. Thus, DEG unigenes had to meet the following criteria: produce a p-value below or equal to 0.05 and a log2-fold value above log2 FC ≥ + 1 for upregulation or produce a p ≤ 0.05 together with a log2fold value below log2 FC ≤ -1 for downregulation.
DEGs were annotated and calculated for enriched ontologies at a significance level of 5% by AgriGO v.2 [80], a specific GO analysis toolkit and a database for agricultural purposes. There are three ontologies in the GO database, namely, molecular function, cellular component and biological process. For this study, we used the results from the category "biological process". DEGs were also annotated using KEGG pathway enrichment analysis aiming at identifying significantly enriched metabolic pathways or signal transduction pathways affected by the NK603 transgene insertion in each variety. Pathways with adjusted p-value < 0.05 were significantly enriched. KEGG analysis was performed using 'clusterprofiler' and 'enrichplot' packages in R environment [86,87]. For the second KEGG analysis an annotation mapping [56,85] was performed using Kobas 3.0 (https ://kobas .cbi.pku. edu.cn/ kobas3). KOBAS is a widely used gene set enrichment analysis tool, its annotation module accepts gene list as input, and generates annotations for each gene based on multiple databases about pathways, diseases, and Gene Ontology.
Selected key genes were examined using DAVID Bioinformatics Resources 6.8 (https ://david .ncifc rf.gov/). After creating a gene list (using Entrez ID), a Gene Functional Classification Analysis was performed. Similarities were measured by Kappa values; furthermore, KEGG pathways were determined. First, a cluster analysis was performed with the following settings: the classification stringency was set to the lowest level, i.e. the individual parameters were as follows: (a) Kappa Similarity: "Similarity Term Overlap level = 3", "Similarity Threshold = 0.2"; (b) Classification parameters: "Initial Group Membership = 3", "Final Group Membership = 3", "Multiple Linkage Threshold = 3". Subsequently, the similarity between all genes was determined using Kappa scores between 0.05 and 1. Cohen's Kappa values are generally accepted to be a robust measure for genetic similarity.
Kappa results range between 0 and 1. The higher the value of Kappa, the stronger the agreement. Kappa more than 0.7 typically indicates a strong agreement between two genes [57].

Confirmation of differentially expressed genes with RT-qPCR
The gene expression of two differentially regulated genes (GRMZM2G127948, AC204711.3_FG003) was assayed by RT-qPCR. To test the biological and the technical reproducibility of the RNA-Seq results (validation of RNA-Seq), a new set of three RNA pools was generated from AG8025conv and AG8025RR2 which was not used for the main RNA-Seq assay. In addition, three new pools of AG9045RR2 and its counterpart AG9045conv were generated to confirm the identified differences between the AG8025conv and AG8025RR2 gene expression.
RT-qPCR was performed on a Rotor Gene Q (Qiagen) using the GoTaq ® 1-Step RT-qPCR System (Promega) according to the manufacturer's instructions. Each reaction was performed in a 20 µL final volume containing 10 ng total RNA and 0.2 µM of each primer. Primers targeting differentially expressed genes were designed using the software Primer-BLAST (https ://www.ncbi.nlm.nih. gov/tools /prime rblas t/). Their sequences are given in the results section. Primer for the reference gene, ubiquitin carrier protein (fwd. 5′-CAG GTG GGG TAT TCT TGG TG-3′, rev. 5′-ATG TTC GGG TGG AAA ACC TT-3′) were taken from Manoli et al. [55]. The primers were purchased from Sigma-Aldrich (Vienna, Austria).
First, to test the specificity of the PCR primer a melting curve analysis was executed which should give just a single peak for each primer pair. Secondly, the efficiency was tested by a dilution series with four different concentrations of a sample. From this test, a standard line and a slope was obtained [68]. Each standard of the dilution was tested in duplicate with a target primer and ubiquitin carrier protein as reference.
All reactions were performed with a hold step at 37 °C (15 min) followed by an initial denaturation step at 95 °C (10 min), followed by 45 cycles of denaturation (95 °C for 10 s), annealing (60 °C for 30 s) and extension (72 °C for 30 s) with a fluorescence measurement at the last step of each cycle. A melting curve, ranging from 60 to 99 °C, with fluorescence measurements at 1 °C intervals, was done after every RT-qPCR, to determine the specificity of the reaction. The differential expression of the two selected genes was measured by comparing the RNA of three independent pools of GM maize grains with three independent pools of conventional maize grains. Each sample pool was tested in triplicate. One pool consisted of 8 maize grains. For inhibition testing and to evaluate the efficiencies of the specific and reference gene PCR assays, standard curves were pipetted in every run. Negative-control and reverse transcriptase (RT)-minus controls (reverse transcription reaction without addition of reverse transcriptase) were also used. Every run was repeated on another day and the mean of the values was taken.
Data were evaluated using the Rotor Gene Q Series software. Linearity (R2) and efficiency (E = 10 [−1/slope] ) of every reaction were within the accepted values. For a valid linearity a value of > 0.98 and for a valid efficiency a slope between −3.1 and −3.6 was required. Relative quantitation was calculated using the 2 −ΔΔCt method [52]. Values were normalized using the reference gene ubiquitin and efficiency correction was performed as described by Pfaffl [64].

1_RNA sequencing (RNA-seq) and library construction for NK603 and near-isogenic maize kernels
Three cDNA libraries per variety were constructed for each GM variety and their near-isogenic non-GM counterparts (12 libraries in total). Table 1 gives an overview of the statistical analyses and the bioinformatic data processing of the libraries per variety. The libraries of AG8025RR2 and AG8025conv yielded 93.6 M of raw reads on average. A total of 14.6% of the raw reads were removed due to adapter removal and cleaning (quality filtering and elimination of rRNA). Next, 73.8% of the raw reads could be mapped uniquely to the B73 maize reference genome version AGPv3. Overall, 5.7% of the raw reads contained repetitive sequences, which aligned at multiple sites. Further, 5.9% of the reads could not be mapped because they had too many errors to match the target sequence. Uniquely mapped reads corresponded to 39,625 unigenes containing 29,146 genes with counts above zero in the transgenic and near-isogenic varieties. The libraries of AG9045RR2 and AG9045conv yielded 71.9 M of raw reads on average; 9.2% of these were removed, 48.6% could be mapped, 4.6% contained repetitive sequences, and 16.4% of the reads could not be mapped. Uniquely mapped reads corresponded to 39,625 unigenes containing 33,420 genes with counts above zero in the transgenic and near-isogenic varieties.

2_Exploratory data-PCA and heatmap
Unsupervised PCA results are depicted in Fig. 2 for the entire data set. The analysis demonstrated a clear cluster by variety (PC1 23% of variation) and a second cluster by transgenic event (PC2 11% of the variation). Variability within the replicates was low (except for one sample from AG9045). For AG9045 and AG9045RR2 the distance between transgenic and non-transgenic is lower than the distance between both non-transgenic. For AG8025 and AG8025RR2 the distance between transgenic and nontransgenic is similar to the distance between both nontransgenic. The heatmap data are presented in Fig. 3. The heatmap hierarchical clustering corroborates the PCA results as the varieties pairs are clustered together. In addition, there was no significant variability between replicates observed.

3_Characterization of significant differentially expressed genes
Based on the criteria outlined in section Materials and methods (subchapter statistical analyses) a total of 2477 genes were determined to be differentially expressed by DESeq2 in AG8025RR2 vs. AG8025conv. (29,146 unigene counts above zero), and a total of 440 genes were determined to be differentially expressed in AG9045RR2 vs. AG9045conv. (33,420 unigene counts above zero). EdgeR analysis indicated 3534 and 694 DEGs were present in AG8025RR2 vs. AG8025conv and in AG9045RR2 vs. AG9045conv, respectively. The DEGs and analytical process results are shown in Fig. 1. The complete list of all DEGs analyzed by both tools are available under Additional files: 1-8. Analysis of the data using EdgeR and DESeq2 showed that a total of 1355 genes were significantly upregulated and 1105 genes were significantly downregulated in the AG8025RR2 variety vs. the AG8025conv counterpart. On the other hand, in the AG9045 variety, a total of 308 genes were significantly upregulated, and 79 genes were significantly downregulated in the AG9045RR2 variety vs. the AG9045conv counterpart. If just common genes were considered, 27 genes were commonly upregulated, and 26 common genes were downregulated in both varieties when comparing the significantly upregulated or downregulated genes of the two varieties AG8025 and AG9045.

4_Singular enrichment analysis, characterization of Gene Ontologies, and REViGO
After statistically defining DEGs, these genes were further analyzed by SEA using the agriGO tool. Exploring GO has become a widespread practice to obtain insights into the potential biological meaning of RNA experiments. For this analysis, there is a growing number of available gene-sets and functional literature containing many genes/proteins and biochemical pathways. Computational analysis helps to find functionally coherent genesets that are statistically overrepresented in a given gene list. The results are mapped as well and compared to gene functional categories of a GO database. Thus, GO analysis may indicate that a certain biological process plays a role in the analyzed biological condition.
For the SEA analysis, the sum of the upregulated plus downregulated DEGs of both varieties, AG8025RR2 vs. AG8025 conv (n = 2460) and AG9045RR2 vs. AG9045conv (n = 387), was calculated (Additional files 9, 10). There were n = 81 GOs in the AG8025RR2 vs. AG8025conv comparison and n = 78-GOs in the AG9045RR2 vs. AG9045conv comparison. However, these numbers were reduced when considering just the   Table 2). Redundancy of GM terms constitutes a major problem for the interpretation of RNAseq results. Very recently, a software called REViGO was developed using an algorithm that evaluates semantic similarities in GO assignments from hierarchical functional annotation of gene ontologies [79]. GO terms with semantic similarities are used to find representative subsets of the terms. Clusters are formed, and each of the GO terms is assigned to the clusters. Cluster representatives are kept, while subordinate cluster members can be removed, thus avoiding redundancy in the results. However, it should be considered that GO annotations are hierarchical, REViGO retains more general ontologies and removes more detailed ontologies, and other programs e.g. ClueGO sort ontologies differently. The usage of REViGO enabled us to reduce the number of ontologies from 24 to 16 GOs. When taking a closer look at the total amount of genes, we discovered n = 81 unique genes (see Table 3) within the 16 ontologies. The results, in particular the number of genes in the GO clusters, the FDR values from the SEA analysis, as well as the p-values for the frequency, uniqueness and dispensability-these parameters are useful to select different level of ontologies-of the SEA analysis are given in Table 2. The largest groups of ontologies were different types of responses, for instance, response to stress, cold, inorganic substances, chemicals, acid chemicals, oxygen-containing compounds, and salt stress, the defense response, and the response to abiotic, biotic, chemical, external or endogenous stimuli. Next, a smaller group of ontologies contained three single-organism processes: a biosynthetic process, singleorganism metabolic process and single-organism cellular process. Other ontologies were developmental processes involved in reproduction and single-organism biosynthetic processes.

5_KEGG pathway classification of identified genes
The GO ontologies evaluated in the previous section contain just a general semantic description of DEGs, which are merged into certain GO terms. One general drawback of GO terms is that their meaning may be very general. However, an important aspect of this study should be that the biological/biochemical functions of identified genes and their functional annotation are understood. This is possible by means of pathway analysis using KEGG annotation/enrichment analysis of single DEGs. KEGG is a database containing a collection of genomes, biological pathways, diseases, drugs, and chemical substances. Thus, KEGG is a popular method to find functionally related genes and pathways that are enriched in a gene list and can be defined based on participation in a metabolic or signaling pathway.
For the KEGG analysis, the two NK603 varieties, AG9045RR2 and AG8025RR2, were compared to their isogenic counterparts, AG9045conv and AG8025conv, respectively. First, an annotation mapping was performed   Table 3). The detailed result of the annotations is specified in Additional file 11. We have also performed a KEGG pathways scatter plot for the enriched metabolic pathways for each of the NK603 varieties. Biosynthesis of secondary metabolites was found altered in both NK603 varieties compared to their near-isogenic conventional counterpart. The variety AG8025RR showed imbalance for starch and sugar metabolism, oxidative phosphorylation, glyoxylate and dicarboxylate metabolism, butanoate, beta-alanine, valine, leucine and isoleucine metabolism and also spliceosome (Fig. 4). Whereas variety AG9045RR only showed arginine and proline metabolism disturbance (Fig. 4). The list of enriched pathways and their corresponding p adjusted values are present in Table 4 and in Additional file 12.

6_DAVID-functional classification of genes
We compared the genes of the overrepresented gene classes (SEA) with the genes involved in significantly affected metabolic pathways (KEGG). Genes occurring in SEA as well as in KEGG analysis of AG8025(RR2) and AG9045(RR2) are depicted in Table 5. If the genes of Table 5 are to be used for more detailed analysis, the knowledge of the genetic similarity between these genes would be very valuable. Thus, the relatedness between the most important genes of Table 5 was determined. Eleven genes occurring in the GO-ontologies as well as in the KEGG pathways of AG8025 and AG9045 were analyzed using the gene functional classification tool of the DAVID Bioinformatics Package (see Materials and methods). The analysis revealed a cluster of three genes (Entrez No. 103631112, 100037774, 542333), while eight genes of the list were not in the output. The genetic similarity between all eleven genes was determined by the calculation of Kappa scores (see Table 6). In addition to the determination of similarities, the three cluster genes were also used to determine Two of these pathways (zma01110 and zma01100) were also found with the above-mentioned KEGG enrichments of all significant DEG genes. Thus, a more accurate representation of the components of these KEGG pathways is intriguing. However, since the KEGG pathways contain a huge number of genes, the number of the background genes occurring in the maize database and the complexity of the maps of KEGG pathways were too high to be listed in the results section, but are indicated in Additional files 14 and 15 instead. Based on this information, it is possible to look up the biochemical relationships between the 11 genes of Table 5.

7_Verification of RNA sequencing results by RT-qPCR
As described above, two different program packages were used for analysis of the primary RNA-seq data with the aim of avoiding false-positive results. However, as in most transcriptomic studies our RNA-seq data contain just a low number of biological replicates (n = 3), but a high number of total data. Thus, from a mathematical perspective it is desirable, and it is practiced as a gold standard to validate DEG data with additional methods such as microarray hybridization or RT-qPCR analysis [18,22,44,49,69,89]. Reverse transcription followed by PCR is a powerful tool for the quantification and detection of gene expression levels, in particular for low-abundance transcripts. Thus, we decided to assign    two randomly selected differentially regulated genes for reverse transcription real-time PCR. One of the genes was upregulated (AC204711.3_FG003), while the other was downregulated (GRMZM2G127948). The results of the DEGs (log2 FC values) tested with DeSeq2 and EdgeR for both genes are depicted in Table 7. The log2 FC values for these genes were either above the cut-off value of log2 FC = + 1 (upregulated) or below log2 FC = −1 (downregulated). The first gene, GRMZM2G127948, codes for the protein ´Caffeoyl-CoA O-methyltransferase 1´. Caffeoyl-CoA omethyltransferase 1 is a key enzyme in lignin biosynthesis. Lignin provides mechanical strength to vascular tissues and protects plants from biotic stresses, including pathogen attack [83]. In investigating the ontologies to which these genes belong, the gene with the ID GRMZM2G127948 could be found in the following ontologies of both varieties-AG9045 as well as AG8025: GO:0044710 (single-organism metabolic process), GO:0044711 (single-organism biosynthetic process), and.
GO:0044763 (single-organism cellular process). The second gene with the ID AC204711.3_FG003 is a gene (senescence-associated protein DIN) relevant for senescence-related processes. Senescence is involved in the deterioration found in several parts of a plant or functional characteristics at a cellular level, in death rates or in fecundity. AC204711.3_FG003 can be assigned consistently to the following functions in the comparison of both varieties-AG9045NK603 vs. AG9045 and AG8025NK603 vs. AG8025: GO:0044710 (single-organism metabolic process), GO:0044711 (single-organism biosynthetic process), GO:0044763 (single-organism cellular process), GO:0009719 (response to endogenous stimuli), GO:0042221 (response to chemicals), GO:0001101 (response to acid chemicals), GO:0006950 (response to stress), and. GO:1901700 (response to oxygen-containing compounds).
To assess if the DGE analysis of both genes can be validated with independent methods, these genes were tested twice using RT-qPCR with the aid of gene-specific and reference primer pairs (see Materials and methods section). The log2 FC data obtained with RT-qPCR compared to the RNA-seq data are shown in Table 7. The log2 FC values of the upregulated gene were distinctly above the positive cutoff value in both GM varieties, log2 FC = + 5.08 and + 2.20 respectively. Additionally, the tendency of the downregulated gene was in accordance with the DGE data. The log2 FC value of the first variety (log2 FC = −1.21) was below the cut-off value (log2 FC = 1), and the log2 FC value of the second variety (log2 FC = −0.82) was very close to the negative cut-off value (log2 FC = −1). For these two genes, the log2 FC values measured by RT-qPCR expression were in good agreement with the DGE data. Thus, the RNA-seq results agreed with the RT-qPCR results and seem to be of high reliability. Nonetheless, it would be still necessary to conduct more single gene experiments using RT-qPCR as a method to get a better estimate of the reliability of the implemented RNA-seq results.

7_Measurement of NK603 grains vs. grains of the near-isogenic conventional variety
When comparing gene expression of GM varieties with those of conventional counterpart plants, it should be considered that potential nonspecific effects might influence the results of the RNA-seq analysis. Thus, we were interested in whether a general difference exists among the investigated kernels. To get a rough estimate of potentially unspecific differences, we measured the weight of GM and conventional kernels. To determine the differences between transgenic and conventional maize, 53 grains of each variety were weighed, and evaluated for significant differences. Because the weights of the grains were not normally distributed (8025conv p = 0.031; 8025RR2 p = 0.00074; Shapiro Wilk test) we evaluated the significance of the differences by a non-parametric test (Mann Whitney U-test). AG8025RR2 maize grains had an average weight of 0.291 g and AG8025conv maize grains 0.378 g. This difference was significantly different (z-value = + 7.12, limit for p = 0.001 significance = + 3.29). Additionally, AG9045RR2 maize grains had an average weight of 0.357 g and AG9045conv maize grains an average weight of 0.232 g. This difference was significantly different as well (z-value = −8.87, limit for p = 0.001 significance = −3.2).
Thus, in the case of the AG8025 comparison, the conventional variety was heavier, whereas in the case of the AG9045 comparison, the transgenic variety was heavier. Figure 5 shows the dry weights of the maize grains of all tested varieties. All measurement values for these samples are present in Additional file 13. As the results were opposite in the two varieties, it cannot be concluded that there is a consistent weight difference between transgenic maize and its conventional counterpart.

Rationale of the study
Gene variation and interactions are common and important phenomena in understanding plant genetics and breeding. Thus, a high-throughput RNA sequencing approach allows a dynamic and functional analysis of maize genetics. The development of 'omics' technology has enabled comprehensive analysis of gene interactions, as well as unintended effects, in different GM events on transcript, protein and metabolite levels [23,27]. However, due to differences in methodological approaches and/or genetic background, little to no consistent results have been obtained among previous studies on how gene expression is influenced by transgene insertion and expression in plant genomes. The rationale of this study was to investigate potential unintended effects deriving from the insertion of a specific transgene-NK603-into two transgenic maize varieties. RNA-seq analysis was used as a molecular profiling technique to study two GM crops in comparison to their near-isogenic maize varieties. Differentially expressed genes were detected with high stringency, taking into account different varieties. False-positives were limited, as we employed different statistical packages for DGE analysis [4,16,53]. Overall, this approach allowed the detection of common effects in two variety pairs. Between the variety pair AG8025RR2 and AG8025conv, we took upregulated DEG genes (n = 2460) that were common between EdgeR and Deseq2 as well as downregulated genes common between EdgeR/Deseq2, of these were 55.1% upregulated and 44.9% downregulated. Between the variety pair AG9045RR2 and AG9045conv, n = 387 genes were differentially expressed (79.6% upregulated and 20.4% downregulated). This number corresponded to ~ 6.21% (for AG8025conv and AG8025RR2) and ~ 0.98% (for AG9045conv and AG9045RR2) of all detected maize genes. Further investigation must be performed to convincingly explain why the variety pair AG8025RR2 has a significantly higher number of differentially expressed genes than the AG9045RR2 variety pair as it does not correlate to the total number of annotated genes in these libraries.

Gene annotations and ontologies
In the analysis of RNA-seq experiments used to study several varieties, it is common to further analyze genes that are affected in both varieties. In our case, it turned out that the 27 upregulated genes and 26 downregulated genes were differentially expressed in the comparison of the AG8025RR2/AG8025conv and AG9045RR2/ AG9045conv variety groups, respectively. To reduce artifacts, we decided to analyze not only two pairs of varieties but also performed several statistical evaluations. For this reason, in addition to individual gene evaluations, we determined matching annotations of GO ontologies and KEGG ontologies for the two GM/ conventional groups. Because each ontology term contains clusters of multiple genes, the ontology results are less affected by single false-positive or false-negative results than DEGs. In the statistical evaluations of the GO ontologies, we used only those ontologies for further analysis that matched between AG8025RR2/ AG8025conv and AG9045RR2/AG9045conv, which consisted of n = 24 ontologies with biological processes in both varieties. A REViGO analysis enabled us to reduce this number to 16 GOs containing n = 81 genes (although other programs such as ClueGO are merging the groups differently). The result of REViGO can be summarized as follows: the largest groups of ontologies were twelve different types of responses. Further ontologies contained three single-organism processes and developmental processes involved in reproduction. When calculating KEGG annotations and even more importantly KEGG enrichments discovering significant pathways with p < 0.05, the results were different between the two variety pairs AG8025(RR2) and AG9045(RR2). Variety pair AG9045(RR2) contained less significant pathways (n = 2) than variety AG8025(RR2) (n = 9). The pathways with p < 0.05 are shown in detail in Additional file 12 and Fig. 4 indicating some main components of the analysis as well as enrichment map of AG8025(RR2) and AG9045(RR2). Several aspects are worth to be emphasized. In variety pair AG9045(RR2) there was a cluster consisting of three pathways (zma00650: Butanoate metabolism; zma00280: Valine, leucine and isoleucine degradation; zma00410: beta-alanine metabolism). It is even more interesting that AG9045(RR2) had a second cluster containing the same pathway found also in AG8025(RR2) (zma01110: Biosynthesis of secondary metabolites).
In order to identify unintended effects arising from the insertion of the NK603 gene more in detail and to facilitate future investigations, we compared the results of SEA/REViGO and KEGG analysis and used the genes that occurred in both evaluations (see Table 5) for additional analyses. Thus, we have determined the similarity among the most important genes with a DAVID analysis and by means of Kappa values. In addition, we have listed the background genes of important KEGG pathways and the biochemical maps of these pathways in Additional files 14 and 15.

RT-qPCR
For a reliable interpretation of RNA-seq analyses, it is necessary to validate the results by additional molecular methods, for example by microarray or by RT-qPCR [3,44,61,69]. La Paz et al. was able to show that a high number of DEGs analyzed with RNA-seq could be confirmed with both microarray and RT-qPCR [49]. We were able to confirm the results of two randomly selected differentially expressed genes using RT-qPCR as well. Evaluation of the analyses shows that the results are highly reliable. These two genes were randomly selected and are important for plant function. One gene with the ID GRMZM2G127948 codes for the protein caffeoyl-CoA O-methyltransferase 1, a key enzyme in lignin biosynthesis. Lignin provides mechanical strength to vascular tissues and protects plants from biotic stresses, including pathogen attack [51,83]. The second gene with the ID AC204711.3_FG003 is a gene (senescence-associated protein DIN) relevant for senescence-related processes. Senescence is involved in deterioration and is found in several parts of a plant or has functional characteristics at the cellular level or may be involved in death rates or fecundity. The expression of these genes can be found in the articles of [16,42,47,59,66,77,78].

Interpretation of genetic diversity using PCA and heatmap
In the interpretation of our results, we are aware that there may be more genes that show differential expression in RNA-seq experiments of genetically modified organisms (GMOs) or when exposed to different environmental conditions. The results of the comparison of single GM crop varieties with corresponding nearisogenic varieties may only apply to the specific variety and to the conditions of the year when the variety was harvested. In order to assess variance structure of our data, we performed PCA and heatmap. The results confirm that we have robust data that differentiates the genotypes even by transgenic event although the highest variability is explained by genotype. AG8025RR2 and its near-isogenic variety have more DEGs than AG9045RR2 and its corresponding near-isogenic variety. Differences between the two near-isogenic varieties (AG9045 and AG8025) are higher than differences between AG9045RR2 and AG9045 but similar compared to the differences between AG8025RR2 and AG8025. This indicates that most of the transgenic effects are variety-dependent which is also implicated by the heatmap results. However, as the two pairs of varieties we examined have different genetic backgrounds, the DEGs occurring in both variety pairs could also indicate that some differences may be conserved among the two varieties. In order to assess the role of GM crops in more detail and to get a better insight into general effects of GM crops, such as unintended effects and pleiotropy [33,50], further studies would have to be performed. For these investigations it would be important to carry out an analysis of additional NK603 varieties containing different genetic backgrounds and the experiments should be performed in more standardized and well defined environmental conditions [8]. To distinguish between environmental and genetic effects, the varieties would have to be cultivated at different times of the year, and the genetic relationship between GM varieties and comparators should be clear in more detail. However, it is difficult to implement this approach, as many varieties, and especially isogenic conventional lines, necessary for RNA-seq analyses of GM crops, are in the hands of corporations and have not been made available for research purposes in response to our enquiries.

Interpretation of unintended genetic effects
Our results indicate that several unintended effects involved in different biological processes have occurred as a result of NK603 integration. Further single-gene analyses of the key genes and genes associated with affected biochemical pathways in different NK603 maize varieties and corresponding counterparts is necessary to interpret the role of the transgene and the biological significance of our findings.
Unintended genetic effects of GM plants have been described in other publications as well [2, 5-8, 12, 17, 19, 32, 35, 40, 45, 49, 58, 60, 84, 88]. According to these studies, transgenic inserts into the genome of one plant variety might affect the overall expression of other endogenous genes in the GM plant. Benevenuto et al. compared proteome profiles of herbicide-tolerant NK603 maize to near-isogenic non-GM maize under drought and herbicide stress and detected twenty differentially abundant proteins mainly assigned to energetic/carbohydrate metabolic processes. When comparing the NK603 maize and its non-GM near-isogenic variety under the same environmental conditions differences were identified in the levels of jasmonate, methyl jasmonate and cinnamic acid and in the abundance of 11 proteins [12]. This is similar to our results as we observed a high abundance of the "jasmonate ZIM domain-containing protein (ko13464)" across the DEGs which, together with their interacting partners, plays an essential role in orchestrating the cross talk between jasmonate and other hormone signaling pathways [89]. Agapito et al. analyzed the proteome profiles of stacked commercial maize hybrid containing insecticidal (BT = Bacillus thuringiensis) and herbicide tolerant traits (NK603) in comparison to the corresponding single event hybrids and non-GM conventional counterparts in the same genetic background as well as in comparison to a non-GM landrace variety under highly controlled growth conditions. Twenty-two proteins were differentially expressed in stacked and single GM events versus non-GM isogenic maize and a landrace variety. These proteins were mainly assigned to energy/carbohydrate and detoxification metabolism. Stacked GM genotypes were clustered together and distant from other genotypes analyzed by PCA. In addition, the varieties containing either BT or NK603 were clustered separately and clearly different from the non-GM varieties [2]. Arruda et al. and Herrera-Agudelo et al. detected significant differences in proteins, metalloproteins, enzymes and metals between transgenic soybeans harboring a RR insert and non-transgenic soybeans [5][6][7]40].
However, many authors of these studies argue that the changes in gene expression, protein distribution and metabolite content caused by genetic modification are less frequent than those caused by environmental factors or natural variations. In the study of Coll et al., natural variation explained most of the variability in gene expression among the samples, but the authors still emphasized that "transcriptional differences of conventionally bred varieties should be considered in the safety assessment of GM plants" [19]. When interpreting these studies, it should be beared in mind that genetic diversity in natural populations can be very large. Breeding or local populations generally show much lower genetic diversity than landraces and different types of wild relatives in maize [34]. Maize has been used as a crop for about 9,000 years [26] and its domestication resulted in a wide genetic diversity of native landraces [82]. During this period, considerable changes in the morphology and physiology of maize may have occurred. Thus, the quantity and quality of genetic effects and the ratio of genetic effects to effects of transgenic insertions may be very diverse in different populations. Specific aspects of risk assessment such as the selection of comparators have been discussed and developed by the EFSA GMO Panel [24].

Secondary effects
Another important requirement for future RNA-seq analyses is to consider secondary effects. La Paz et al. described a small but significant delay in seed and plant maturation, which possibly influenced the functional annotation and expression of differentially expressed genes of MON810 plants [49]. To consider secondary effects in our analysis, we measured the weight of maize grains and compared them between NK603 and isogenic varieties. In these experiments, we did not find any consistent effects between the two varieties. However, we do not know whether there are other secondary effects between GM crops and conventional varieties. Thus, it would be important to carry out additional experiments, for which one could either evaluate certain parts of plants or otherwise analyze embryos as performed in the study of La Paz et al. [49].

Conclusion
Overall, our data clearly demonstrate substantial differences between the analyzed transgenic varieties and their non-transgenic counterparts. PCA confirms a distinct difference between conventional and transgenic samples for variety AG8025 and a slight difference for the other variety pair. Several biological processes and metabolic pathways are modulated in the transgenic varieties. These differences indicate that several unintended effects have occurred as a result of NK603 integration. Heatmap data imply that most of the transgenic insert effects are variety-dependent. However, identified key genes involved in affected pathways of both variety pairs show that transgenic independent effects cannot be excluded. Further research is necessary to clarify the role of internal and external influences on gene expression.
In general, our studies show that transcriptomic analysis is very useful to assess gene interactions, pleiotropic effects and unintended effects in transgenic crops. Thus, this technique may be a valuable tool for assessing genes that affect plant health or plant fitness, and serve as complementary safety analysis for the pre-market approval of GMOs. Even though RNA-seq has become the standard method for transcriptome analysis there are still analytical gaps that need to be taken into account, especially those related to the quantification of low levels transcripts [22]. However, with continuous developments of RNA-seq strategies, it is anticipated that more robust transcript identification will be able to be performed from longer reads. Thus, allowing a more accurate detection of individual, allele-specific biological variations and splice variants [61,81].
In future investigations of NK603 varieties, it would be important to analyze genes of affected biochemical pathways more in detail to assess the influence of internal (such as variety or NK603 transgene) and external (such as environment) factors on the expression of these genes. Moreover, ontology clusters of this study should be evaluated more accurately and with other methods for example, different types of omics studies [48,58]. Furthermore, analyzing stacked varieties [84] and the consideration of several transgenic events could be of great interest for the regulatory authorities.