Systematic miRNome profiling reveals differential microRNAs in transgenic maize metabolism

Background While some genetically modified organisms (GMOs) are created to produce new double-stranded RNA molecules (dsRNA), in others, such molecules may occur as an unintended effect of the genetic engineering process. Furthermore, GMOs might produce naturally occurring dsRNA molecules in higher or lower quantities than its non-transgenic counterpart. This study is the first to use high-throughput technology to characterize the miRNome of commercialized GM maize events and to investigate potential alterations in miRNA regulatory networks. Results Thirteen different conserved miRNAs were found to be dys-regulated in GM samples. The insecticide Bt GM variety had the most distinct miRNome. These miRNAs target a range of endogenous transcripts, such as transcription factors and nucleic acid binding domains, which play key molecular functions in basic genetic regulation. In addition, we have identified 20 potential novel miRNAs with target transcripts involved in lipid metabolism in maize. isomiRs were also found in 96 conserved miRNAs sequences, as well as potential transgenic miRNA sequences, which both can be a source of potential off-target effects in the plant genome. We have also provided information on technical limitations and when to carry on additional in vivo experimental testing. Conclusions These findings do not reveal hazards per se but show that robust and reproducible miRNA profiling technique can strengthen the assessment of risk by detecting any new intended and unintended dsRNA molecules, regardless of the outcome, at any stage of GMO development. Electronic supplementary material The online version of this article (10.1186/s12302-018-0168-7) contains supplementary material, which is available to authorized users.


Background
Genetically modified (GM) crops intended for market release have been mostly created through in vitro DNA modification to encode a novel protein. However, a new generation of genetically modified organisms is being designed to change their RNA content in order to regulate gene expression [1]. The regulation of gene expression by RNA molecules, specifically double-stranded RNA (dsRNA), is now known to be possible through a process called RNA interference (RNAi). The biochemical pathway for this type of gene regulation can be generally understood as the disruption of the production of proteins when there are dsRNA molecules interacting with the DNA sequence of that gene. RNAi is an important biological pathway that is used by many different organisms to down-regulate the expression of endogenous or exogenous mRNA by the interaction of small interfering RNAs (siRNAs) and microRNAs (miRNAs) with catalytic RNA-induced silencing complex (RISC) to bind to homologous mRNA sequences and to cleave them [2]. miRNAs are known to control a wide range of biological processes in eukaryotic organisms, from embryogenesis [3], to developmental and cellular behaviour [4], and response to biotic [5] and abiotic stress [6]. miRNAs function as gene regulators and they have been shown to play important roles in multiple plant Open Access *Correspondence: sarah.agapito@genok.no † Sarah Zanon Agapito-Tenfen and Vinicius Vilperte contributed equally to this work and their names follow alphabetical order 1 GenØk-Centre for Biosafety, Forskningsparken i Breivika, Sykehusveien 23, 9294 Tromsø, Norway Full list of author information is available at the end of the article developmental and signalling pathways through small RNA biogenesis or related pathways [7]. Due to the environmental stability of dsRNAs, a property perhaps overlooked based on the relative instability of single-stranded species of RNA [8], GM crops can be engineered to produce dsRNAs that are harmful to insects and worms that feed on these plants [8][9][10][11]. On the other hand, studies have claimed that dsRNA molecules may rapidly degrade and dissipate in aquatic microcosms [12], and are biological inactive in agricultural soils after 2 days of exposure [13]. Nonetheless, the stability and transmissibility of dsRNAs suggest the potential for the existence of exposure routes that are relevant to environmental and human risk assessments of GM organisms [14]. In such cases, all new intended and unintended dsRNA molecules should be identified in the GM organism and products [15]. This should be obtained through a semi-targeted qualitative profiling of small RNA molecules, using next-generation sequencing in a comparative assessment of the GM and near-isogenic non-transgenic variety [14]. Since dsRNA produced by plants can target important mRNAs from pathogens feeding in the crop and cause lethality, there is a growing interest in using RNAi for pest and disease control, both as a traditionally applied insecticide/fungicide and within GM plants [16]. However, while some GMOs are intended to produce new regulatory-RNA molecules, these may also arise in other GMOs not intended to express them [14,15,17]. GM crops are regulated by international and national legislation and many countries adopt pre-market risk assessment to evaluate any risks that GM plants may pose to animal and human health and the environment [18,19]. There are similarities and differences in the risks associated with RNAi-based GM plants relative to those posed by GM crops producing insecticide and herbicidetolerant proteins [14]. As part of the pre-market risk assessment, many regulatory authorities evaluate all relevant scientific data on the molecular characterisation of the GM plant in question, such as the source and function of the donor DNA, the transformation method, the organization of the inserted DNA at the insertion site(s), and the expression and stability of the insert [18,19].
Our previous paper documented risk assessment advice offered to government regulators during official risk evaluations of GM plants for use as human food or for release into the environment, how the regulator considered those risks, and what that experience teaches us about the GMO risk assessment framework [14]. We concluded that the process should include inter alia experimental procedures that would identify all new intended and unintended dsRNA molecules in the GM product and we have suggested improvements for such.
The European Food and Safety Authority opinion is in agreement with our recommendations [1].
Since the evidence that exogenous plant miRNAs are present in the sera and tissues of various animals and that these exogenous plant miRNAs are primarily acquired orally through food intake [20], the safety of GMOderived miRNAs has been debated [21][22][23][24][25]. Nonetheless, further studies have detected the presence of exogenous miRNAs in human samples from uncommon dietary sources (e.g. rodents), indicating that exogenous miRNAs originate from technical artefacts rather than dietary intake [26]. Hence, there are no conclusive evidences and consensus about this topic.
Only few efforts designed to investigate unique risk assessment issues related to GM crops or products containing dsRNA have been reported. The present study is the first that characterizes the miRNome of commercialized GM maize events to further investigate potential alterations in miRNA regulatory network. We were interested in testing miRNA content in two of the most commonly used transgenes in agriculture which currently do not produce intentional dsRNA; the insecticide Bt and the herbicide-tolerant Roundup Ready transgenes. The seed set of stacked (combination of transgenes by crossing or mating) and single GM maize events, as well as the non-transgenic near-isogenic counterpart, developed in the near-isogenic genetic background, enables the isolation of potential effects derived from transgenes independently.
Here, we show that the expression profile of several miRNAs is altered in GM varieties and that novel and potential transgenic miRNA sequences were detected. Although no unequivocal explanation was found for these phenomena, the interactions between miRNA, their targets and biochemical pathways are discussed. These analyses were carried out in order to create a documented record of the emergence of such alterations. The experimental approach of this study may also expand the knowledge about the potential risks of RNAi-based GM plants through the identification of potential changes that require further characterisation of products from cutting-edge technology.

Deep sequencing of sRNA species in transgenic maize varieties
In order to identify conserved maize miRNAs in transgenic maize varieties, sRNA libraries were constructed from leaf material and sequenced using Illumina highthroughput technology. We have performed miRNA profiling in a unique set of stacked and single transgenic maize hybrid varieties, as well as the conventional counterpart, which were all developed under the same genetic background. Each sample consisted of a pool of 10 plants and three biological replicates were analysed, with a total of 30 plants assessed per variety. A flowchart of the miRNA analysis pipeline is available in Fig. 1. The full sequencing dataset is available through Sequence Read Archive (SRA) at NCBI under identification number SRP160077.
After removing low-quality sequences, the total library size ranged from 13,306,004 to 19,618,997 reads, with 44-54% belonging to the 18-26 nt range and thus used for further analyses. Reads with sizes < 18 and > 26 nt were excluded from the datasets. The read size distribution for all libraries is available in Additional file 1. The majority of the sRNAs with 18-26 nt belong to 24 nt size category, accounting for 13.2-15.9% of the total number of reads, followed by 21 nt (6.1-8.0%) and 22 nt (4.7-6.0%) (Fig. 2a). Similar distribution patterns were also observed in previous studies investigating plant miRNAs [27][28][29][30].
On average, approximately 22% of the reads matched to other types of non-coding sRNAs, such as rRNAs, tRNAs, snRNAs or snoRNA, and approximately 18% matched organellar RNA (chloroplastidial and mitochondrial) ( Table 1). In addition, we have built a reference sequence database with the transgenic insert sequences of MON89034 and NK603 events, which was also used   for mapping the reads. Transgenic sequences were detected in all GM varieties. The stacked variety showed a higher number of sequence reads (4908 reads) which was close to the sum of the reads in both single event varieties (476 read for the Bt variety and 4173 reads for the RR variety). Further investigation of the potential for those sequences to be an active miRNA was performed using miRPrefer prediction tool. The analysis did not yield significant results because the sequences did not fulfil the criteria for either precursor sequence size and/ or maximum of 2 nt 3′ overhangs in mature-star duplex.
No reads matching the transgenic sequences were found in control samples. The B73 RefGen_v3 database was used as reference genome to search for conserved miRNA sequences. The overall distribution of miRNAs is represented by a Venn diagram in Fig. 2b. A total of 128 miRNAs were found to be shared among Bt, RR, RRxBt and control sample libraries, accounting for ~ 90% of the annotated miRNA detected in our bioinformatics analysis. Other miR-NAs (10%) were assigned to be present in one variety only, these are: zma-miR2118c in control, zma-miR395f in RRxBt, zma-miR395g in Bt, and both zma-miR395d and zma-miR395 h in RR samples. Nonetheless, these exclusive miRNAs were only present in a few numbers of reads. Our sequencing effort presented a good miR-Nome coverage with a total of 404,728 to 428,487 number of reads per variety matching to 143 known miRNAs sequences from 27 miRNA families in maize (out of known 156 miRNA present on the reference genome) (Fig. 2c). The largest family was zma-MIR169 with 17 members, followed by zma-MIR171 (13 members) and zma-MIR156 and zma-MIR395 (11 members).
In relation to the abundance of each miRNA family, zma-MIR166 was the most present in all varieties (50.5 to 55.2% of the total miRNA content), followed by miR168 (14.5-16.2%) and miR156 (7.8-8.8%). Out of the top ten expressed miRNA families, our study shared 6 miRNA families matching the miRNA maize reference paper [30]. A heatmap with the normalized read count for all the identified miRNA families for each variety is shown in Fig. 2d. Although single event varieties group separately from stacked and control samples in the heatmap, there was no statistically distinct pattern of miRNA families' abundance between varieties. The list of all identified maize miRNAs and the respective normalized read count for each variety is available in Additional file 2.
Defining transgenic maize-specific microRNA signatures through multiple microRNA profiling comparisons Pairwise differential expression analysis was performed based on the normalized read count using DESeq2 R package for each identified miRNA. To facilitate visualization, we have used the logarithm base 2 of the fold-change (Log2FC) as a measure of accumulation levels. Therefore, a gene up-regulated by a factor of 1.5 has a Log2FC of 0.58, a gene down-regulated by a factor of 1.5 has a Log2FC of − 0.58, and a gene expressed at a constant level (with a FC of 1) has a Log2FC equal to zero [31].
Thirteen known maize miRNAs were found to have statistically significant differences in abundance (p-adjusted < 0.05 was used as the only analytical threshold) ( Table 2) between maize varieties. The results of the differentially abundant miRNAs were validated using SYBR reverse transcriptase quantitative PCR (RT-qPCR). The amplification of miRNA sequences is based on forward primers designed to match the full mature miRNA sequences, whereas the reverse primer is universal and targets poly-A tail added to miRNA sequences (Additional file 3).
Noteworthy, four miRNAs belonging to the same family could not be distinguished due to their identical mature miRNA sequence. We also observed that the amplification of the zma-MIR399eij showed a considerable difference when compared to the RNA-Seq results. This might be due to the homology of the designed primer to other members of the family, such as zma-MIR399c and zma-MIR399d, which might have hampered the results. In addition, zma-miR162 and zma-miR529 showed multiple melting curve peaks for the RT-qPCR analysis, thus precluding a reliable quantification of their expression. Overall, the comparative cycle threshold levels generated for most of the selected miRNAs were the same as those determined by Illumina sequencing, indicating that the sequencing data produced in this study were accurate. The comparison of Illumina and RT-qPCR quantification results is shown in Table 3.
One of the challenges in elucidating the biological functions of miRNAs is to identify their regulatory targets. We have predicted targets for the 13 differentially regulated miRNAs using psRNAtarget and found 116 potential endogenous transcript targets. These included splicing variants and gene expression inhibition by either mRNA cleavage or translational inhibition. The list of all miRNA targets is available in Additional file 4. The miR-NAs with the largest number of mRNA targets were zma-miR529-5p (47) and zma-miR169f-5p (25).
As expected, many of the putative targets of maize conserved miRNA identified were transcriptions factors (TFs) mRNAs, which is considered a general trend for plant miRNAs [7]. Of the 116 predicted targets, 24 are members of transcription factor gene families (i.e. member of the nuclear transcription factor Y subunit A-3 family, squamosa promoter-binding protein-like (SBP domain) transcription factor family, putative GATA transcription factor family and the Ocs element-binding factor 1-putative bZIP transcription factor superfamily) involved in many central developmental and physiological processes, including morphogenesis, abiotic and biotic stress responses.
The miRNome profile of Bt-expressing transgenic plants was the most distinct one showing nine out of 13 dys-regulated miRNAs from five miRNA families (zma-MIR162, zma-MIR167, zma-MIR169, zma-MIR399 and zma-MIRNA529). We have further focused our downstream analysis of the Bt miRNome to explore possible factors influencing miRNA accumulation in this variety. We were able to verify a statistically significant increase of cry1A.105 protein transcript accumulation in Bt versus stacked event through RT-qPCR (Fig. 3). Although these results do not directly correlate to the most altered miRNA profile in the Bt variety, the expression and accumulation of CRY proteins suggest alterations in the miR-Nome of both single and stacked varieties, which showed the two most distinct miRNA profiles. However, because CRY protein in a non-native maize protein, it is not possible to verify annotated pathway networks that integrate CRY.
In order to gain insight into what metabolic pathways would be affected by the altered miRNAs in the Bt variety, we submitted transcript targets of these differentially regulated miRNAs to a singular pathways enrichment analysis using agriGO toolkit. We observed that these transcripts were mainly assigned to RNA biosynthetic process (GO:0032774), regulation of DNA-dependent transcription (GO:0006355) and regulation of RNA metabolic process (GO:0051252) (Fig. 4a). Examples of mRNA targets associated with these pathways are dicer-like 1 protein family, Pumilio-family RNA binding domain, Histone H4.3, Nucleic acid binding protein NABP, Nuclear transcription factor Y subunit A-3 isoform 2, squamosa promoter-binding protein-like (SBP domain) transcription factor family protein, putative GATA transcription factor family protein, among others. These proteins shared common molecular functions of transcription regulator activity (GO:0030528, GO:0003700) and nucleic acid binding (GO:0003676, GO:0003677) ( Fig. 4b and Table 4).
Stacked transgenic event plants had the second most distinct miRNome with four miRNAs (zma-miR399e, zma-miR399i, zma-miR399j and zma-miR827) being differentially accumulated when compared to control or to single event expressing CP4-EPSPS. The three zma-MIR399 species found to be less accumulated in stacked event target the same mRNA assigned to contain an HSF-type DNA-binding domain. On the other hand, zma-miR827 targets 13 different gene transcripts, among them transcripts containing domains of WD40-like Beta Propeller Repeat, BTB/POZ domain, SPX-Major Facilitator Superfamily domain and Myb/SANT-like DNA-binding domain (Additional file 4).
We have performed independent gene target profiling using RNA-Seq analysis to evaluate the expression of predicted target transcripts in all varieties and the non-transgenic control (Additional file 5). The RNA-Seq differential analyses showed 77, out of 116, putative miRNA-predicted target relations. Although there is no established threshold for the miRNA-target ratio, computational prediction and experimental validation have demonstrated negative correlations among miRNAtarget pairs due to mRNA degradation which were also observed in this study [32].
Overall, the results showed that 13 miRNAs were differentially modulated in different pairwise comparisons. The results also suggest that CRY proteins can be considered a major factor influencing the miRNome of both Bt and stacked events in an additive manner. The altered miRNAs can target a range of endogenous transcripts, mostly transcription factors and nucleic acid binding domains that are key molecular functions for basic genetic regulation and some have been confirmed by independent RNA-Seq analysis.

Prediction of novel miRNA candidates and their potential targets
In order to identify novel miRNAs sequences in our samples, we have used the miR-PREFeR pipeline, together with MiPred online tool, to predict their pre-miRNA structure. A total of 20 putative novel miRNAs have been identified in one or more varieties (Table 5). Detailed information of novel miRNAs is available in Additional file 6. The novel miRNAs sequences were temporarily named as 'zma-MIR-number' followed by an 'X' to differentiate them from the conserved maize miRNAs (e.g. zma-MIRX01). The predicted structure of all 20 novel miRNAs is available in the Additional file 7, and three of these novel miRNAs are shown in Fig. 5. Out of the 20 novel miRNAs, two were exclusively detected in control samples (zma-MIRX03 and zma-MIRX14), one in the RR samples (zma-MIRX11) and one in the Bt samples (zma-MIRX13).
Confirmatory analysis of novel miRNA was performed by RT-qPCR for those which was possible to design primers (Additional file 3). We were able to obtain single amplicons by melting curve analysis for only 7 miRNA primers. All these 7 miRNAs (zma-miRX03, zma-miRX05, zma-miRX06, zma-miRX10, zma-miRX11, zma-miRX13 and zma-miRX20) have provided positive amplifications via RT-qPCR (Table 5). Surprisingly, amplifications were also obtained in samples that showed no miRNA reads by miRNome sequencing. Spurious amplifications might be explained by the amplification of similar or identical miRNA mature sequences originated from different pri-miRNA loci. Poor sequencing of the sample and/or limitations of the bioinformatics analysis is also a possibility. In fact, in spite of the availability of a few amplificationbased platforms for miRNA detection, the common technical challenges still exist, namely sensitivity, specificity especially for one base variation, selectivity between precursor and mature miRNAs [33].
The prediction of the novel miRNA targets showed 420 potential endogenous transcript targets, including splicing variants and gene expression inhibition by either mRNA cleavage or translational inhibition. The list of all novel miRNA targets is available in Additional file 8. Transcripts targets of miRNAs shared by all samples were submitted to singular pathways enrichment analysis using agriGO toolkit. Only two GO terms (p value < 0.05) were significantly enriched for biological processes in the miRNA targets of the predicted novel miRNAs, cellular lipid metabolic process (GO:0044255) and lipid biosynthetic process (GO:0008610). Four major molecular functions had significant GO term hits and most transcripts were assigned to catalytic activity (GO:0003824) ( Fig. 6 and Additional file 9). The two novel miRNAs identified in control samples only (zma-MIRX03 and zma-MIRX14) had identical targets and those were mainly identified as Det1 complexing ubiquitin ligase (DDA1), putative calciumdependent lipid-binding (CaLB domain) family protein

Identification of miRNA variants or "isomiRs"
The miR-PREFeR pipeline was used to analyse possible variant sequences arising from the miRNA precursor. We considered as variants, or isomiRs, any sequence that mapped to the miRNA precursor and was neither identified as being a mature nor a star sequence by the pipeline. Out of 156 annotated miRNAs in maize genome (B73, RefGen_v3, release 25), the algorithm was able to predict 96 miRNA sequences, with all of them showing at least one variant sequence. In Fig. 7, we present two cases of such miRNA variants. In the first example, the corresponding miRNA* is more abundant than the annotated miRNA (Fig. 7a). In the second example, the most abundant sRNA is neither the annotated miRNA nor the miRNA*, but one of its variants (Fig. 7b).

Discussion
In this report, we aimed at identifying which correlation of the differentially expressed miRNAs is more likely to play roles in transgenic maize through miRNA-gene interaction networks compared to non-transgenic nearisogenic control. Unless the GM plant is intended to produce new dsRNA molecules, the RNA itself is rarely formally considered in a risk assessment [14]. This is surprising since many GMO risk assessment guidance draw special attention to the characterization of novel RNAs, frequently mentioning the need to provide information on any expressed substances in the recombinant-DNA plant, such as: "any gene product(s) (e.g. a protein or an untranslated RNA)" (paragraph 32 of [34]); "demonstrate whether the inserted/modified sequence results in intended changes at the protein, RNA and/or metabolite levels" (page 10 of [35]); "detailed description of the expression of the gene product inserted in the host organism" (item 10, Annex II [36]). Although the GM varieties studied here were not tested for their dsRNA content during approval process because of assumptions made either explicitly or implicitly in the context of risk assessments, our detailed analysis of MON89034, NK603 and stacked MON89034xNK603 suggests that their miRNA profiles are not equivalent to the near-isogenic non-transgenic control samples.
We have developed a study design helping to fill the vacuum for the risk assessment of dsRNAs unique to, or present at specific concentrations in, GMOs. We have also identified areas of research that need technical improvements as well as areas of current uncertainty   related to specific aspects of miRNA biology that require further investigation. Our study may be useful in the first step of a sequential approach to assess the potential for adverse effects arising from dsRNA-initiated modifications to organisms before bioinformatics and exposure analysis are performed. We have sampled GM plants  grown under strictly controlled conditions and we have compared their miRNA profile to unmodified nearisogenic control lines that have the same genetic background as the GM varieties so to isolate the differences in miRNA content only. Profiling approaches have been long recommended by EFSA for the allergenicity assessment of the whole GM plant [37]. In addition, a recent report from EFSA's workshop on 'Risk assessment considerations for RNAi-based GM plants' also highlighted the need for more research on metabolic profiling of such plants to support risk assessment through the use of various 'omics' techniques [1]. Although comparative miRNA profiling has not yet been performed in commercialized GM crops, other layers of genetic information, such as the proteome and metabolome, have been widely applied and shown to be a useful source of information for the risk assessment. Many studies have investigated the proteomes, transcriptomes and metabolomes in grains and leaves of insect-resistant and herbicide-resistant maize and suggested metabolic alterations due to the genetic transformation process [38][39][40][41][42]. Due to a more recent wave of omics datasets that have started to be used in some risk assessment areas and have been accepted as a powerful tool to substitute or complement current studies, EFSA has also hosted a dedicated colloquium early this year to move towards the implementation of omics techniques to GMO risk assessment in Europe [43]. Remarkable findings of ingested plant miRNA in animal liver and blood have been demonstrated by a Chinese study and this directed the attention of GMO regulators to potential cross-kingdom regulation between GM plants and higher organisms. Plant microR-NAs have been detected in human blood, demonstrating that dsRNAs can survive digestion and be taken up from the gastrointestinal tract [20]. The plant-derived dsRNA molecules silenced an endogenous gene in human tissue culture cells, and in mouse liver, small intestine and lung [20]. Other studies have also shown that regulatory RNAs from plants can be detected in animals, including humans [24,[44][45][46], and it was found that some dsRNAs Fig. 7 Distribution of sRNAs along conserved miRNA precursors. a Example where the corresponding miRNA* is more abundant than the annotated miRNA; b example where the most abundant sRNA is not the annotated miRNA or the miRNA*, but one of its variants from plants were more frequently present in animal tissues and blood than predicted from their level of expression in plants. Other plant-derived dsRNAs (i.e. rice, corn, barley, tomato, soybean, wheat, cabbage, grapes and carrot) have later on been detected in humans [47]. Although the functionality of such exogenous RNAs in terms of mediating gene expression in animals has been controversial [48,49], there is enough evidence that such RNAs are capable of translocating from the host to its interacting organism, and vice versa [50]. Overall, there is an enormous lack of knowledge about the biological significance of the uptake of plant-derived miRNAs and, therefore, further investigations are required [49].
In the current study, we performed miRNA expression profiles, as well as miRNA target prediction analysis to identify a subset of miRNAs that could be correlated with transgenic maize varieties that expressed insecticidal CRY protein and the EPSPS protein that confers herbicide tolerance. The transgenic plants encoding two CRY toxins, event MON89034, revealed the most distinct miRNome profile followed by the profile of the stacked transgenic event. Yet, the distinct miRNome is not correlated to a higher expression of CRY transcripts in Bt single event variety. Our independent RT-qPCR analysis revealed similar expression of Cry2Ab2 and a significantly decreased expression of Cry1A.105 in single compared to the stacked variety. In fact, quantification of transgenic transcripts and its correlation to protein accumulation has frequently yielded inconsistent results [51].
We have then examined the relationship of altered miRNAs to pathways that are perturbed by the corresponding miRNA gene targets in Bt samples. Nine miR-NAs from five miRNA gene families were exclusively up-regulated in Bt plants (except for zma-MIR399 which was down-regulated in Bt samples). Metabolic pathway enrichment analysis revealed that relevant gene targets were related to RNA metabolic processes and regulation of DNA-dependent transcription. These observations suggest built-in redundancy of overlapping gene targets regulated by different miRNAs targeting the same set of pathways. The potential for functional redundancy among plant miRNAs has been previously studied by inactivation of the entire MIR164 family in Arabidopsis [52]. It was found that mir164abc triple-mutant plants showed severe defects in flower development and phyllotaxis, and these defects were not observed in plants in which only individual MIR164 genes are disrupted. Thus, there were indications that miR164 miRNAs are controlling shoot development in a redundant manner, while the degree to which individual miR164 miRNAs contribute to the regulation of different developmental processes varies [52]. Therefore, differences in the expression patterns of the individual miRNA genes imply that redundancy among them is not complete and that these miRNAs show functional specialization, which may pose extra challenges for any single miRNA to be used as a target for further studies on its potential biosafety implications.
Another altered metabolism in Bt samples was related to primary metabolic process (GO:0080090 and GO:0044238) and nitrogen compound metabolism (GO:0051171 and GO:0006807). Primary metabolism in plants is directly involved in growth, development, and reproduction, involving basic processes such as carbohydrate, lipid and protein metabolism. Our previous proteomic investigation using the same GM varieties also showed differential modulation of enzymes that catalyse chemical reactions involved in carbohydrate metabolism [40]. Overall, other comparative proteomic studies of transgenic versus non-transgenic crops corroborate these results. In fact, the energetic metabolism, including the carbohydrate metabolism, has been the most frequently observed metabolism category within comparative proteomic analysis of transgenic versus non-transgenic crops (see compilation at Table 3

from [41, 53]).
Two miRNAs families had several miRNAs being altered, zma-miRNA167 in Bt plants and zma-miRNA399 mostly in stacked plants. zma-miR167 is known to target auxin response factor (ARF) genes [54]. The phytohormone auxin influences many aspects of plant development and the identity of several miRNA targets suggests roles for miRNAs in auxin signalling [55]. ARFs are a plant-specific family of DNA binding proteins that bind to auxin-responsive promoter elements (AuxREs), which are found in early auxin response genes, including Auxin/ Indole-3-Acetic Acid (Aux/IAA), small auxin-up RNA (SAUR), and auxin-inducible GH3 genes, and can either enhance or repress transcription [56]. The role of auxin dys-regulation in Bt plants is yet unknown but might be related to the CRY protein pathway as part of the plant defence system against pathogens [57]. zma-MIR399 was less abundant in stacked plants and showed a target transcript that contains a plant calmodulin-binding domain. Calmodulin-binding proteins are thought to be involved in diverse pathways that are dependent on Ca 2+ signalling, such as plant development and adaptation to environmental stimuli [58]. In addition, also less abundant, zma-MIR827 targets another class of transcription factors, such as the bZIP transcription factor superfamily, the Ocs element-binding factor 1 and the Rossmann-fold superfamily protein. In plants, basic region/leucine zipper motif (bZIP) transcription factors regulate several processes including pathogen defence and stress signalling [59].
Off-target effects of ncRNAs producing plants have been observed. Tomato plants expressing short-hairpin RNAs targeting a replicase protein of Tomato leaf curl New Delhi virus (ToLCNDV) showed phenotypic abnormalities, such as needle-shaped leaves, reduced stature, and decreased lateral rooting, among others [60]. Similarly, soybean plants containing RNAi construct aimed to silence soybean myo-inositol-1-phosphate (GmMIPS1) gene have also shown off-target effects, such as impaired development of seeds [61]. A study conducted by the product owner assessed the NK603 event as 'substantially equivalent' to its isogenic counterpart by a nutrient composition analysis (e.g. amino acids, fatty acids, vitamins and minerals) [62]. On the other hand, when 'omics' techniques were used to profile the proteome and metabolome of the same event, alterations in energy metabolism and oxidative stress pathways were found [53]. miRNA profiling methodologies are still evolving and responding to new information and hence, no consensus framework or pipelines have been established. We have chosen to use a more comprehensive analysis applying RNA-Seq as high-throughput unbiased method of miRNA target identification to uncover networks that require large-scale analysis instead of demonstrating individual miRNA:mRNA interactions (i.e. by RT-qPCR) because it misses the capacity for miRNAs to regulate complex gene networks [63]. We were able to confirm 77 out of 116 putative miRNA-predicted targets by this approach. Inconsistent correlation might be related to data observed amongst a pool of indirect changes in transcript abundance and, although it may assist in describing the predominant genes and pathways affected by a miRNA, it does not distinguish between direct targets [63]. It could also be related to the different target prediction algorithms, which are based around a model of interaction between the 5′-end of the miRNA called the 'seed region' and the 3′ untranslated region (3′-UTR) of the mRNA and this might not be the only interaction ruling between a miRNA sequence and its target. Indeed, increasing evidence demonstrates that targeting can also be mediated through sites other than the 3′-UTR and that seed region base pairing is not always required [63]. Nonetheless, identifying direct targets remains problematic given the potential modest effect on levels of some target transcripts and the fact that some miRNA targeting might occur predominantly at the level of translational repression [64].
Twenty novel miRNA sequences were found in the different maize samples. Two of them were observed in control samples only, one in RR sample only and another one in Bt samples alone. The other 15 novel miRNAs were found in all samples in similar abundance. The identification of novel miRNAs is a common approach in many plant miRNome profiling studies due to the yet relatively small-annotated miRNA database for plants [65]. Novel miRNA confirmatory analyses revealed amplifications for 7 miRNA primers but no differential profiles were observed. Target prediction of novel miRNA is possible through the use of bioinformatics database search according to sequence homology with known genes. Our target prediction revealed approximately 420 gene transcripts that could be regulated by these miRNAs. Metabolic pathway enrichment analysis of the 15 novel miRNAs shared among all samples showed GO terms that were assigned to lipid metabolism. Although several bioinformatics tools and algorithms have been developed to identify novel miR-NAs in silico, there is still a need for detailed in vivo analyses to confirm their existence and expression. This is also the case for transgene-derived 22-24nt RNA sequences. We have found several potential transgenic miRNA which could not be verified using in silico prediction tools. Sequence variation in cis-regulatory elements which regulate genome transcription by binding to transcription factors might lead to higher expression of the associated miRNA locus only in the GM varieties. Whether such variation exists and is caused by transgene itself or the transgene-insertion process or other unrelated mutation, it needs additional experiment.
In addition to the novel miRNA sequences and the potential transgenic miRNA sequences, we found isomiRs for 96 miRNAs in our dataset. Frequently overlooked, isomiRs are several length and/or sequence variants of the same annotated miRNA. These variants were originally dismissed as experimental artefacts [66]. Because they vary in RNA length and/or sequence, many isomiRs are capable of showing major effects on miRNA targeting efficiency, AGO incorporation, loading into the RISC complex and stability [66].
We have previously alerted that unanticipated offtarget adverse effects can be difficult to detect and it is not yet possible to reliably predict them using bioinformatics techniques [14]. Therefore, we aimed at identifying all unintended dsRNA molecules in the GM product through a semi-targeted qualitative profiling of small RNA molecules using next-generation sequencing in a comparative assessment between the GM and conventional parent. This approach would fit perfectly into the first phase of the molecular characterization step of current GMO RA. Nonetheless, identified unintended changes in miRNA content should be further verified and tested according to the guidelines we have previously suggested for improvements in risk assessment of GM crops or products containing dsRNA [14].
In this study, we have used a profiling approach which was able to identify several unintended miRNA sequences in single and stacked GM maize. These unintended miRNA sequences arose from differential expression of conserved miRNA families with a significant impact on metabolic pathways, novel miRNA sequences not yet annotated in the maize database, miRNA variants of unknown function (isomiRs) and potential transgenic miRNA sequences. Due to technical limitations, some of these sequences could not be verified or tested in silico and/or in vitro and therefore, require further in vivo analysis. All of these newly identified sequences can be a source of potential off-target effects and should be further considered in risk assessment frameworks [14]. In addition, safety studies cannot only rely on bioinformatics to predict potential miRNA gene targets. Our study showed that many of the predicted targets had no annotation in the maize GenBank and thus testing would be needed to specifically assess potential adverse effects on animal and human health and the environment.

Conclusions
In summary, our data revealed 13 miRNAs that may play important roles in transgenic maize plants through the dys-regulation of gene targets. These 13 miRNAs generally modulate unique gene targets but co-ordinately target a few common pathways, with partially overlapping key nodal molecules, highlighting the redundancy of miRNAs and gene targets that focuses on reprogramming critical pathways to facilitate transgene expression in maize cells. Thus, the complex gene networks regulated by these miRNAs highlighted the challenges in validating these results and also to guiding hypothesis-driven downstream analysis. Twenty novel miRNAs, miRNA variants (isomiRs), and potential transgenic miRNA sequences were also detected in our experiment.
We have generated data from a level of information that has never been considered before in risk assessment of genetically modified plants reaching the market. Production of intended dsRNA molecules in GM crops is a legitimate concern due to their potential off-target effects on silencing genes other than those intended. We have provided an approach to properly assess the safety of dsRNA-producing GM organisms through the identification of all new intended and unintended dsRNA molecules in the GM product by profiling omics technique. Nonetheless, we have pointed to technical limitations and to areas in which further experimental validation is required. Experimental procedures will no doubt evolve as new information becomes available, but this should not be a reason not to apply our current scientific and technological knowledge to assessing the safe use of dsRNA in a commercial scale.

Plant material and growth chamber conditions
We have chosen commercial transgenic maize varieties tolerant to the herbicide glyphosate ("Roundup Ready", RR) and Bt-expressing plants that have been extensively cultivated in Brazil and other southern hemisphere countries like Argentina and South Africa over the last years. Four maize varieties were used in this study. All are commonly found in the Brazilian market: AG8025RR2 (unique identifier MON-ØØ6Ø3-6 from Monsanto Company, glyphosate herbicide tolerance, Sementes Agroceres), AG8025PRO (unique identifier MON-89Ø34-3 from Monsanto Company, resistance to lepidopteran species, Sementes Agroceres) and AG8025PRO2 (unique identifier MON-89Ø34-3 × MON-ØØ6Ø3-6 from Monsanto Company, stacked event resistant to lepidopteran species and glyphosate-based herbicides, Sementes Agroceres). As a control, we have used the non-GM near-isogenic variety AG8025 (Sementes Agroceres). AG8025 seeds are hybrid progeny of the single cross between maternal endogamous lines "A" with the paternal endogamous line "B". In this case, all hybrid progeny is 100% AB genotype. The stacked variety MON-89Ø34-3 × MON-ØØ6Ø3-6 expresses two insecticidal proteins (CRY1A.105 and CRY2AB2 proteins derived from Bacillus thuringiensis (Bt), which are active against certain lepidopteran insect species) and two identical EPSPS proteins, under the control of distinct regulatory sequences, providing tolerance to glyphosate-based herbicides. More information on the generation of these commercial transgenic varieties and their expression cassettes can be freely obtained at Biosafety Clearing House hosted by the UN Convention on Biological Diversity (http://bch.cbd.int/), searching by the variety unique identifiers provided above.
Seedlings were grown side by side in growth chambers (EletrolabTM model 202/3) set to 16-h light period and 25 °C (± 2 °C), using Plantmax HT substrate (Buschle & Lepper S.A.) and watered daily. No pesticide or fertilizer was applied. Around 50 plants were grown in climate chambers out of which 30 plants were randomly sampled per maize variety. The collected samples were separated into three groups of ten plants. The ten plants of each group were pooled and were considered one biological replicate. Maize leaves were collected at V4 stage (20 days after seedling). Leaf pieces were cut out, weighed and placed in 3.8-ml cryogenic tubes before immersion in liquid nitrogen. The samples were kept at − 80 °C until RNA isolation.

miRNA isolation and deep sequencing
Total RNA was extracted from approximately 100 mg of frozen leaf tissue using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The isolated RNA was quantified using Nan-oDrop 1000 (Thermo Fisher Scientific, Wilmington, USA) and resolved in 1% MOPS denaturing gel. RNA samples (1 μg) were sent to FASTERIS SA (Geneva, Switzerland) for library construction and sequencing. Twelve small RNA (sRNA) libraries were constructed. Before library construction, samples were quantified with Qubit ® Fluorometric Quantitation (Life Technologies, California, USA) and quality control was accessed with Agilent 2100 Bioanalyzer (Agilent Technologies, California, USA). Briefly, the construction of sRNA libraries consisted of the following successive steps: (i) acrylamide gel purification of RNA fragments corresponding to sRNAs sizes (18-30 nt); (ii) ligation of the 3p and 5p adapters and indexes to the RNA in two separate subsequent steps, each followed by acrylamide gel purification; (iii) synthesis of cDNA followed by another acrylamide gel purification; and (vi) final step of PCR amplification to generate a cDNA colony template library for Illumina sequencing.

Library analysis of small RNAs
All low-quality reads with FASTq values below 13 were removed, and 5′ and 3′ adapters, as well as index sequences, were trimmed using the Genome Analyzer Pipeline (Illumina) at Fasteris SA. In addition, quality control was performed with FastQC tool (http://www. bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) and no bases with Phred quality score below 30 were found; thus no additional trimming was necessary. Reads outside the 18-26 nt range were excluded from the analysis. sRNAs belonging to rRNAs, tRNAs, snRNAs and snoRNAs, as well as chloroplastidial and mitochondrial sequences, derived from Zea mays and deposited in the Ensembl Plants and NCBI GenBank databases, were identified through mapping using the Bowtie2 v.2.2.4 software [67] and removed from the dataset.

Identification of conserved and novel miRNAs
In order to determine conserved maize miRNA, the filtered sRNA sequences were mapped to the Zea mays genome (B73, RefGen_v3, release 25) deposited in the Ensembl Plants database using the Bowtie2 v.2.2.4 software [67], in which the complete alignment of the sequences was required, and no mismatches were allowed. In addition, reads with multiple alignments on the genome were further excluded from the analysis. The prediction of novel miRNA was performed using the miR-PREFeR pipeline [68] with the following parameters: (1) Maximum length of 250 nt for a miRNA precursor; (2) Reads depth cut-off of 200; (3) Maximum gap length of 50 nt between two contigs to form a candidate region; (4) Minimum and maximum length of the mature sequence of 21 and 24 nt, respectively; (5) No requirement of the star sequence to be expressed; and (6) Allow the mature star duplex to have only 2 nt 3′ overhangs. All annotated mRNA and miRNA in the B73, RefGen_v3, release 25 were excluded from the analysis. In addition, each putative novel miRNA was searched against the maize miR-NAs deposited at miRbase (http://www.mirba se.org/). Moreover, the web-based tool MiPred [69] was used to classify the predicted novel miRNA as real, pseudo or not a miRNA precursor. To be considered a real precursor, the miRNA sequence has to have a Minimum Fold Energy (MFE) < − 20 kcal/mol and a p-value < 0.05. As for the pseudo precursor, at least one of the conditions has to be true.
The miR-PREFeR pipeline was also used to analyse possible miRNA variant sequences for the conserved miR-NAs. The following parameters were used: (1) Maximum length of 300 nt for a miRNA precursor; (2) Reads depth cut-off of 2; (3) Maximum gap length of 100 nt between two contigs to form a candidate region; (4) Minimum and maximum length of the mature sequence of 18 and 24 nt, respectively; (5) No requirement of the star sequence to be expressed; and (6) Allow the mature star duplex to have only 2 nt 3′ overhangs. An annotation file containing 156 conserved maize miRNAs was provided in order to analyse only previously annotated miRNAs.

Differential expression of miRNA
Read counts were retrieved using SeqMonk v0.29.0 (http://www.bioin forma tics.babra ham.ac.uk), in which probes were designed around known miRNA sequences in the reference genome. Only the reads that exactly overlapped the probes were considered for the read counting. The statistical analyses were performed through R language and environment [70] using the DESeq2 R package v1.8.1 [71] to normalize the read counts and perform the differential expression analysis. DESeq2 uses the Wald test for pairwise comparisons, as well as Benjamini-Hochberg method [72] for the p-value adjustment for multiplicity of testing.

Validation of miRNAs by RT-qPCR
Validation of miRNA sequencing was performed by RT-qPCR amplification using designed primers targeting all differentially expressed miRNAs and novel miRNAs found in this study (Additional file 3). Poly(A) Tailing and cDNA synthesis was performed using the NCode ™ VILO ™ miRNA cDNA Synthesis Kit (Invitrogen, California, USA) according to the manufacturer's recommendations. Forward primers consisted of the entire sequence for each of the mature miRNAs, while only one universal reverse primer, provided with the commercial kit, was used. The universal reverse primer annealed in the Poly(A) tail of the cDNA, ensuring that only Poly(A)-cDNA was quantified. Endogenous reference genes were based on previous work on validation [73]. The primers used for the endogenous genes are the same ones used in our previous work, but only the forward primer. The amplification efficiency was obtained from relative standard curves provided for each primer and calculated according to Pfaffl equations [74]. For each of the biological replicates, three independent cDNA syntheses were performed. The RT-qPCRs were performed using the EXPRESS SYBR ® GreenER ™ miRNA qRT-PCR kit (Invitrogen) according to the manufacturer's recommendations. The normalized relative quantity (NRQ) was calculated according to the Pfaffl equations [74]. Information on real-time data for this study has followed guidelines from the Minimum Information for Publication of Quantitative Real-Time PCR Experiments [75].

Relative quantification of transgene transcripts
Quantification of cry1A.105, cry2Ab2 and espsps transcripts were performed on the same samples used for miRNA and mRNA sequencing. Reverse-transcription quantitative PCR (RT-qPCR) assay was adapted from previously developed assays for the specific detection of MON-89Ø34-3 × MON-ØØ6Ø3-6 transgenes [76] to hydrolysis ZEN-Iowa Black ® Fluorescent Quencher (ZEN/IBFQ) probe chemistry (Integrated DNA Technologies, INC Iowa, USA). cDNA was synthesized and amplification of each target gene was performed using the QuantiTect Probe RT-PCR Kit (Qiagen) according to the manufacturer's instructions. The methodology (e.g. PCR conditions, primer and probe sequences, statistical analysis) was the same as described for our previous work [40].

Prediction of miRNA targets and pathway enrichment analysis
The prediction of gene targets for the differentially expressed miRNAs was performed using the psRNAtarget online tool [77]. The following parameters were applied: (1) Maximum expectation: 3.0; (2) Length for complementarity scoring (hspsize): 20; (3) Target accessibility-allowed maximum energy to unpair the target site (UPE): 25; (4) Flanking length around target site for target accessibility analysis: 17 bp upstream and 13 bp downstream; and (5) Range of central mismatch leading to translational inhibition: 9-11 nt.

Confirmatory analysis for miRNA-predicted target relation
Total RNA was also isolated from the same samples using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions and submitted to transcriptome sequencing at FASTERIS SA (Geneva, Switzerland). Sequencing of the cDNA libraries was conducted using the HiSeq SBS Kit v4 (Illumina ® ) in an Illumina HiSeq 2500, with number of cycles of 2 × 125 + 7 (paired-end) in one lane of the HiSeq Flow Cell v4 (Illumina ® ). Base-calling was performed using the pipelines HiSeq Control Software 2.2.38, RTA 1.18.61.0, CASAVA-1.8.2.
All low-quality reads with FASTq values below 13 were removed, and 5′ and 3′ adapter, as well as index sequences, were trimmed using the Genome Analyzer Pipeline (Illumina) at Fasteris SA. Moreover, quality control with FastQC tool (http://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) was performed and no bases with Phred quality score below 30 were found; thus no additional trimming was necessary.
The filtered reads were mapped against the Zea mays genome (B73, RefGen_v3, release 25), deposited in the Ensembl plants database [80], using the TopHat2 v2.1.0 tool [81], with a maximum of two mismatches, in which gaps count as mismatches. Pairwise differential analyses were conducted using the Tuxedo package [82], which uses Cufflinks v.2.2.0, Cuffmerge and Cuffdiff tools, and is able to assemble transcripts, estimate their abundance, and test the expression and differential regulation of RNA-Seq libraries.