GmDREB1 overexpression affects the expression of microRNAs in GM wheat seeds

MicroRNAs (miRNAs) are small regulators of gene expression that act on many different molecular and biochemical processes in eukaryotes. To date, miRNAs have not been considered in the current evaluation system for GM crops. In this study, small RNAs from the dry seeds of a GM wheat line overexpressing GmDREB1 and non-GM wheat cultivars were investigated using deep sequencing technology and bioinformatic approaches. As a result, 23 differentially expressed miRNAs in dry seeds were identified and confirmed between GM wheat and a non-GM acceptor. Notably, more differentially expressed tae-miRNAs between non-GM wheat varieties were found, indicating that the degree of variance between non-GM cultivars was considerably higher than that induced by the transgenic event. Most of the target genes of these differentially expressed miRNAs between GM wheat and a non-GM acceptor were associated with abiotic stress, in accordance with the product concept of GM wheat in improving drought and salt tolerance. Our data provided useful information and insights into the evaluation of miRNA expression in edible GM crops.


Introduction
The global hectarage of genetically modified (GM) crops has increased 100-fold, from 1.7 million hectares in 1996 to 179.7 million hectares in 2015, for 17 to 18 million farmers in 28 countries.The latest data for 1996 to 2014 shows that GM crops contribute to food security, sustainability and climate change [1].However, in many areas, cultivation of the GM crops was prevented because of the dominant ideological voices of the opponents, debates on their potential negative environmental impact, and their health risks to consumers.Although there is no sufficient scientific experimental evidence [2], the adverse effects of several GM crops that have been published in scientific journals still generate suspicion and controversy in the opinion of both the public and the scientific community [3][4][5].Therefore, research assessing the potential risks associated with GM crops is essential for fast scientific adoption and general public acceptance.The current evaluation system for GM crops is focused on proteins, fats, carbohydrates, toxins, and nutritional ingredients.However, miRNAs have not yet been taken into account in this system of evaluation.It was recently found that several gma-miRNAs are differentially expressed in GM soybean seeds compared with the non-GM control [6].miRNAs have been widely demonstrated to play fundamental roles in gene regulation in most eukaryotes [7][8][9][10].The first cross-kingdom regulation of gene expression between humans and edible plants by miRNAs was reported in 2012 [11].Moreover, it was recently demonstrated that Brassica miR-NAs can regulate the expression of human genes and proteins in vitro [12], which has raised new concerns regarding the miRNA components of GM crops.
miRNAs regulate gene expression by guiding target mRNA cleavage or by translational inhibition [7][8][9][10].The targets of miRNAs include a large set of transcription factors (TFs) [13], such as Auxin Response Factors (ARFs) and QUAMOSA Promoter Binding Protein-Like (SPLs), which are involved in the regulation of miRNA expression in a negative feedback loop in Arabidopsis [14].On the other hand, miRNAs are processed from hairpin precursors (pre-miRNAs), and pre-miRNAs are derived from primary transcripts (pri-miRNAs) transcribed from genomic DNA which is most likely regulated by TFs.An over-expressed transcription factor (TaDREB3) in barley affects the expression of miRNAs and other small non-coding RNAs in its leaf [15].
Due to their involvement in the regulation of many stress-related genes, many TFs were utilized in GM crops to improve crop tolerance to abiotic stresses.GM drought tolerance is an extremely important goal given that droughts will likely become more severe and more frequent as climate change impacts crop productivity, agriculture and society.So far, only 3 GM drought tolerant plants have been recommended for commercialization approval including GM maize in America, soybean in Argentina, and sugarcane in Indonesia.GM DroughtGard™ tolerant maize, first planted in the US in 2013, increased more than 15-fold from 50,000 hectares in 2013 to 275,000 hectares in 2014 and 810,000 hectares in 2015, reflecting a high farmer acceptance, at a 3-fold, year-to-year increase between 2014 and 2015 [1].Wheat (Triticum aestivum L.) is one of the most important cereal crops worldwide, and it is used as a stable food grain and as a primary source of straw for animal feeding.Wheat production is severely affected by various abiotic stresses such as drought, salinity, low temperature, and heat, which results in an estimated 50-60% loss in grain yield annually [16].Thus far, except for the glyphosate herbicide tolerance wheat (Mon71800), no other GM drought tolerant wheat has been approved for commercialization.Transgenic wheat transformed with a DREB TF (GmDREB1) from soybean (Glycine max L.), which is driven by an ubiquitin promoter showed strong tolerance to drought and salt [17,18].
In this study, we investigated small RNAs from the dry seeds of a GM wheat line overexpressing GmDREB1 and the non-transgenic control varieties using deep sequencing technology and bioinformatic approaches.Comparative analysis showed that 7 known tae-miRNAs were differentially expressed in the dry seeds between GM wheat and the non-GM acceptor.Combined with miRcheck prediction and experimental verification, 16 conserved and novel tae-miRNA candidates were also identified.This indicates that overexpression of the transcript factor GmDREB1 affects the expression of miRNAs in wheat seeds.Most of the target genes of the differentially expressed miRNAs between GM wheat and the non-GM acceptor were associated with abiotic stress, in accordance with the product concept of the GM wheat line T349 [17,18].These results provided useful data for a bio-safety assessment of GM crops and valuable information for wheat miRNA research.

GM wheat samples and RNA isolation
The transgenic wheat line T349, non-transgenic acceptor Jimai 19 (J19) and comparable nontransgenic varieties Jimai 22 (J22) and Lumai 21 (L21) were provided by the Institute of Crop Science, Chinese Academy of Agricultural Sciences.The transgenic wheat line T349, transformed with GmDREB1 from soybean, showed strong tolerance to drought and salt stresses [17,18].The transgenic and non-transgenic wheat varieties were planted in a field with watersaving measures in the Jinan province of China.The mature wheat seeds were harvested for RNA isolation.The main wheat cultivars, J22 and L21, were cultivated in the same climatic and geographic ecological region as J19.The total RNA was isolated from wheat seeds that were pooled from 3 individuals of the same transgenic or nontransgenic line using the Plant RNA Extraction kit (TaKaRa MineBEST) according to the manufacturer's instructions.The quantity and quality of total RNA were evaluated using a NanoDrop 2000 Spectrophotometer and an Agilent 2100 RNA 6000 Nano kit.

Small RNA library construction and sequencing
The sRNAs were isolated from the total RNA using the TruSeq Small RNA Sample Preparation kit.The sRNA was ligated with a 3'-linker and electrophoresed on a 15% polyacrylamide gel.The 15-30 nt sRNAs were excised from the gel and recovered in 0.3 M NaCl.The recovered sRNAs were then ligated with a 5'-linker and subject to reverse transcription PCR.The RT-PCR products were separated on a 3.5% agarose gel.The 140-160 bp products were selected for library construction and submitted for sequencing on a Hi-Seq 2500 analyzer.The raw data from the small RNA libraries were deposited in the NCBI Sequence Read Archive (SRA) under accession No. SRP091415.

Bioinformatic analysis of high-throughput data
For deep sequencing data analysis, similar methods reported in a previous work were used in this section [6].In brief, raw reads were cleaned by quality and adaptor trimming.Identical reads were collapsed and recorded for abundance.The SOAP programme was used for genome mapping.Only perfectly matched sRNAs with a length between 15-30 nt were subjected to further analysis [19].The wheat (Triticum aestivum L.) genome was downloaded from the cerealsDB database (http://www.cerealsdb.uk.net/cerealgenomics/CerealsDB/copyright.php).

Novel miRNA prediction
Before novel miRNA prediction, small RNA reads were searched against the Rfam database (V11.0) to remove the known non-coding RNAs [22].Customized Perl scripts, RNAfold and miRCheck programmes were used for context extraction, secondary structure determination and novel miRNA prediction [23].Novel miRNA candidates were clustered with known miR-NAs, including mature and precursor sequences, to find novel members in conserved families.Details of this process are described in a previous study [24].

Target prediction and function analysis
Target genes were predicted using psRNATarget [25], a plant sRNA target analysis server (http:// plantgrn.noble.org/psRNATarget).Eight degradome libraries (GSM1606478, GSM1606479, GSM911923, GSM911924, SRR1197125, SRR1197126, SRR1197127, SRR1197128) were downloaded and analyzed for validating the cleavage sites in target genes using CleaveLand4 [26].To further understand the function and classification of the predicted miRNA targets, GO classification of the target genes was conducted with WEGO web service (http://www.geneontology.org/),GO terms assigned to the query sequences and catalogued groups were produced based on their biological process, molecular functions, and cellular components.

RT-PCR assays
Differential expression of the tae-miRNAs was verified by qRT-PCR experiments.One microgram of total RNA was reverse-transcribed and amplified using miRcute miRNA First-Strand cDNA synthesis kit (Tiangen, KR201) and the miRcute miRNA qPCR detection kit (Tiangen, FP401).The comparative ΔΔCT method was used for the relative quantitation of these tae-miRNAs [27].miR159 was selected as the endogenous reference gene for normalization [28].For the validation of novel tae-miRNA candidates, a diluted cDNA template was used for RT-PCR assays with specific forward primers and a universal reverse primer.The amplified products were detected by electrophoresis on a 2% (w/v) agarose gel and validated by Sanger sequencing.
The expression profiles of some target genes of the differentially expressed miRNAs were analyzed by qRT-PCR analysis of a randomly selected set of genes.DNaseI-treated RNA was used for first-strand cDNA synthesis using Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and oligo (dT)18 primers according to the manufacturer's protocol.The comparative ΔΔCT method was used for the relative quantitation of genes [27].The Tubulin gene was amplified as an internal control.The primers used in this section are listed in S1 Table.

Deep sequencing of small RNAs in wheat seeds
Twelve small RNA libraries from T349, J19, J22 and L21 wheat cultivars, each with 3 replicates, generated 187,919,197 raw reads via the Illumina/Solexa deep sequencing platform.After the removal of low-quality reads and 3' adaptor trimming, 159,723,105 (100%) clean reads, corresponding to 25,123,347 (100%) unique reads, were obtained for subsequent analysis.Genomic mapping results showed that 53,820,203 (33.7%) clean reads and 11,937,579 (47.5%) unique reads could find at least one perfectly matching locus in the wheat genome.We regarded 40,599,438 (25.4%) clean reads and 11,669,683 (46.4%) unique reads as unknown reads after searching against the miRBase and Rfam databases, which were used for novel miRNA prediction (Table 1).The length distribution of the clean reads showed that both 21-nt and 24-nt classes were dominant in the total reads (Fig 1A) and only the 24-nt class occupied the majority in the unique reads with obvious superiority (Fig 1B).

Expression profile of known tae-miRNAs in transgenic and nontransgenic wheat
To detect the known tae-miRNAs, the clean reads were aligned against the known tae-miRNAs represented in miRBase 21.As a result, 7 tae-miRNAs were up-regulated in the seeds of GM wheat T349 compared with the non-GM acceptor J19.Moreover, no down-regulated miRNAs were detected (Table 2).The log2 ratios varied from 1.01 to 1.20 and the maximum value was 1.20 for tae-miR319.Among these 7 known tae-miRNAs, 3 miRNAs, tae-miR167c-5p, tae-miR156 and tae-miR9661-3p were up-regulated only in GM wheat seeds compared to a non-GM acceptor, while there were no significant differences between the different non-GM wheat varieties tested in this study (L21/J19 and J22/J19).Four other tae-miRNAs were also differentially expressed between the non-transgenic wheat seeds.Compared with J19, tae-miR319 was up-regulated both in T349 and L21; tae-miR164 and tae-miR9678-3p were upregulated both in T349 and J22 and tae-miR9659-3p was up-regulated in T349, L21 and J22 (Table 2).Meanwhile, the expressions of these known tae-miRNAs were confirmed by qRT-PCR (Fig 2).The differentially expressed miRNAs between the non-transgenic wheat varieties were not listed here, if they were not differentially expressed between T349 and J19.

Predicted novel miRNAs in GM and non-GM wheat
In the 12 libraries, the sequenced miRNAs that were not reported in Triticum aestivum but were in Aegilops tauschii, Brachypodium distachyon, Oryza sativan, Sorghum bicolor or Zea mays according to miRBase 21, were defined as conserved miRNAs.The sequenced miRNAs that were not reported in Triticum aestivum or the other five monocotyledons referenced above were defined as novel miRNAs.Sixteen conserved miRNAs and 49 novel miRNAs were found in wheat seeds after sequence alignment (Fig Compared with the non-GM acceptor J19, 5 conserved miRNAs and 11 novel miRNAs were differentially expressed in the seeds of GM wheat T349 (Table 3).One novel miRNA was  down-regulated and another 15 miRNAs were up-regulated in T349.In addition to T349, the conserved miRNAs ata-miR396c-5p, osa-miR164c and bdi-miR827-3p were also up-regulated in J22.The miRNAs ata-miR396c-5p and ata-miR166a-3p were also up-regulated in L21 and J22 compared with J19 (Table 3).Among the 11 novel miRNAs families, 2 miRNAs were differentially expressed only in GM wheat seeds, while no significant differences were observed between the different non-transgenic wheat varieties tested in this study (L21/J19 and J22/J19).
We defined these two novel miRNAs as tae-N1 (5'GCCTCCGTAGCATAGTGGT3') and tae-N2 (5'GCGTCTGTAGTCCAACGGT3'). The other 9 miRNA families were differentially expressed in both GM wheat T349 and the non-GM varieties J22 and/or L21 compared with J19.The expressions of some novel and conserved miRNAs were tested by qRT-PCR ( 2.4 Target genes of up-and down-regulated miRNAs in GM and non-GM wheat miRNAs regulate target genes involved in plant development and stress response.miR167 target auxin response factor mediated plant root development [29,30].miR156 targets SQUA-MOSA promoter binding protein-like (SPL) TFs related to flowering time, phase change and leaf initiation rate [31,32].miR319 targets TCP transcription factors and miR164 targets NAC domain gene, which contributes to leaf development [33].Moreover, these target genes are all involved in stress response [34][35][36][37][38] (Table 4).To reveal the function of the predicted novel miRNAs and the conserved miRNAs, we predicted their targets using psRNATarget (http:// plantgrn.noble.org/psRNATarget/).As is common for wheat, some predicted targets were functionally unknown and most miRNAs were found to target more than one gene (S3 Table ).Part of these cleave sites were validated by wheat degradome data (S4 Table ).T-plots of the identified targets were shown in S1 Fig.
The expression profiles of some target genes of differentially expressed miRNAs were analyzed by qRT-PCR, and some target genes were down-regulated in T349 compared with J19 (S2 Fig) .GO analysis results of all the target genes of the differentially expressed miRNAs in GM and non-GM wheat showed that many targets were related to water transport and responses to oxygen-containing compounds (S3 Fig) , suggesting that the corresponding miRNAs may be involved in drought/salt stress and oxidative stress.

Over-expressed transcription factor (GmDREB1) in wheat affects the expression of miRNA in seeds
The global area for the cultivation of transgenic plants has grown continuously during the past two decades, despite constant controversy concerning their unforeseeable biosafety issues [1].
The application of transgenic technology has led to significant changes in modern agriculture and this trend seems irreversible [1].Because of the dynamic property of RNA molecules, it is reasonable to propose that some endogenous RNAs might be diverse in transgenic plants.To date, very few reports have discussed the RNA ingredient content difference between GM and non-GM plants [6,15].Herein, we investigated the differential expression of miRNAs in wheat seeds from a GM line (T349), its non-GM receptor (J19) and local major wheat cultivars (J22 and L21) using high throughput sequencing technology.As a result, GmDREB1overexpression in wheat affects the expression of miRNAs in seeds.Other reports also mentioned the relationship between miRNAs and DREB.The expression levels of DREB1A and DREB2A were relatively induced in miR408 overexpressed plants compared to the vector control upon drought stress [39].An over-expressed transcription factor (TaDREB3) in barley affects the expression of miRNAs and other small non-coding RNAs in its leaf [15].In this study, the degree of variance, including miRNAs (Fig 2) and sRNA variance (S4 Fig), between cultivars was much higher than that induced by the transgenic event.Twentythree differentially expressed miRNAs were found and confirmed between the T349 and J19 lines.Moreover, some of the miRNAs affected by DNA manipulation overlapped with the differentially expressed miRNAs related to varietal influence.Further analysis showed that most target genes of these differentially expressed miRNAs were associated with abiotic stress.This is understandable because a GmDREB1 transcription factor was introduced in the transgenic line T349 to elevate its drought and salt tolerance ability [17,18].https://doi.org/10.1371/journal.pone.0175924.t004

Target genes of the differentially expressed miRNAs were associated with abiotic stress
Differentially expressed miRNAs were observed between the GM wheat T349 and the non-GM acceptor J19.Several known miRNAs (miR319, miR164, miR167, and miR156) are involved in drought response.miR319 regulates TCP transcription factors [33] and the NAC domain gene, which contributes to leaf development and drought stress response [38].Moreover, the NAC domain gene is known as the target gene of miR164 in rice [40], and TCP genes have been shown to regulate miR164 in Arabidopsis [34].It was reported that miR319 also target MYB transcription factors [41] that play significant roles in stress responses.miR167 regulates its target gene, ARF, which mediates plant root development [30,37] and responds to drought stress [35].Another identified target of miR167, IAA-Ala Resistant3 (IAR3) has the same roles as ARF in drought stress and root development [29].miR156 targets SQUAMOSA promoter binding protein-like (SPL) TFs, which control flowering time, phase change and leaf initiation rate [31,32,42].The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants.Blocking the miR156 signalling pathway with 35S::MIM156 (via target mimicry) increased the sensitivity of the plant to stress treatment, whereas the overexpression of miR156 increased stress tolerance [36].In transgenic barley overexpressing the drought tolerant gene TaDREB3, miR156 was over two-fold up-regulated when comparing with the non-transgenic control.Thus these known miRNAs up-regulated in GM wheat may enhance the stress tolerance of transgenic wheat.In another study, miRNA microarray analysis showed that miR156, miR167, miR164, miR319, miR396 and miR166 were up-regulated in leaf or root of bread wheat under drought stress [43].GO analysis of all the target genes of differentially expressed miRNAs in the GM and non-GM wheat lines show that many targets are related to water transport and response to oxygen-containing compounds (S3 Fig) , suggesting that the corresponding miRNAs may be involved in drought/ salt stress and oxidative stress.

Possible environmental risks of differentially expressed miRNAs in GM wheat
A controversial and attention-drawing issue is whether the exogenous miRNAs of GM crops could regulate human gene expression by cross-kingdom RNAi, which has been observed in both animal and plant systems.In animal systems, there are a few instances of cross-kingdom gene silencing.One example is that RNAi can be induced by soaking Caenorhabditis elegans worms in RNA solutions or by feeding them antisense RNA-expressing bacteria, such as Escherichia coli [44].Plants transfer RNAi signals into interacting organisms, such as filamentous fungi, oomycetes, nematodes, parasitic plants, and pests [45][46][47][48] to silence their genes; however, it is not confirmed whether miRNA mediates human-plant interactions.Although it was reported that exogenous rice MIR168a is present in the human and mouse sera and the exogenous MIR168a inhibits human/mouse LDLRAP1 gene expression in liver [11], many subsequent reports have also questioned this controversial result [49][50][51][52].It was recently demonstrated that Brassica miRNAs could regulate the expression of human genes and proteins, but only in vitro [12].In this study, the targets of 23 differentially expressed miRNAs in GM wheat seeds were predicted in human (S5 Table ) and chicken (S6 Table ), both of which may eat GM wheat seeds.The top 50 predicted targets with the highest scores of 23 differentially expressed miR-NAs in GM wheat seeds were mainly involved in the regulation of metabolic process, including primary, macromolecule and RNA metabolic processes in humans (S5A Fig) Moreover, they were primarily involved in epithelial cell differentiation in chicken (S5B Fig) .However, our study did not conclude whether there are possible environmental risks of eating GM wheat seeds exist until it was the cross-kingdom RNAi was proven to be induced between animal/ human-plant in vivo, and the differentially expressed miRNAs target these genes as predicted.
In addition to miRNA, sRNA includes other classes, such as small interfering RNAs, a second class of small regulatory RNAs that can direct DNA methylation and regulate salt tolerance and disease resistance [53,54].Moreover, long noncoding RNAs that emerged as important regulators of many biological processes in animals were identified in higher plants [55].In this study, besides the 23 miRNAs discussed above, other sRNAs may be taken into account during GM crop evaluations in the future.The function of these different types of sRNA will gradually be identified in future studies.Moreover, the targets of some miRNAs were not predicted.For example, tae-N1 and tae-N2 from this study may function by interacting with other mRNA or lincRNA.The complications of sRNA justify taking sRNA into account in the evaluation system for GM crops.Our findings provide useful data for assessing the potential risks associated with GM crops and provide valuable information for wheat miRNA research.

Fig 2 .Fig 3 .
Fig 2. Verification of differentially expressed tae-miRNAs by qRT-PCR.The comparative ΔΔCT method was used for the qRT-PCR experiments and miR159 was selected as the reference.https://doi.org/10.1371/journal.pone.0175924.g002 Fig 2).The 11 novel miRNAs and the 5 conserved tae-miRNA candidates were subjected to experimental verification by RT-PCR and Sanger sequencing, and they were all positively identified in at least one library (Fig 4, S3 File).

Table 4 .
Target prediction of known tae-miRNAs in wheat seeds.