- Research article
- Open Access
Pesticide application has little influence on coding and non-coding gene expressions in rice
BMC Genomics volume 20, Article number: 1009 (2019)
Agricultural insects are one of the major threats to crop yield. It is a known fact that pesticide application is an extensive approach to eliminate insect pests, and has severe adverse effects on environment and ecosystem; however, there is lack of knowledge whether it could influence the physiology and metabolic processes in plants.
Here, we systemically analyzed the transcriptomic changes in rice after a spray of two commercial pesticides, Abamectin (ABM) and Thiamethoxam (TXM). We found only a limited number of genes (0.91%) and (1.24%) were altered by ABM and TXM respectively, indicating that these pesticides cannot dramatically affect the performance of rice. Nevertheless, we characterized 1140 Differentially Expressed Genes (DEGs) interacting with 105 long non-coding RNAs (lncRNAs) that can be impacted by the two pesticides, suggesting their certain involvement in response to farm chemicals. Moreover, we detected 274 alternative splicing (AS) alterations accompanied by host genes expressions, elucidating a potential role of AS in control of gene transcription during insecticide spraying. Finally, we identified 488 transposons that were significantly changed with pesticides treatment, leading to a variation in adjacent coding or non-coding transcripts.
Altogether, our results provide valuable insights into pest management through appropriate timing and balanced mixture, these pesticides have no harmful effects on crop physiology over sustainable application of field drugs.
Insect pests (IPs) are the most prominent threats in achieving global food demands of a rapidly growing population. IPs affect the latent yield of all agricultural crops either directly or indirectly. Direct damage may include deformations or necrosis of plant tissues or organs, fouling and dispersion of plant pathogens, while the loss of harvest quality and increase in the cost of crop production may involve indirect damage [1, 2].
Owing to the severity of agricultural insects’ problem, it has become a great challenge to use sustainable measures to control IPs that could affect crop yield. Various control strategies including mechanical, biological, cultural, transgenic and chemical have been followed by farmers to manage IPs since past. Modern biotechnology and genetic engineering led to the development of Genetically Modified Organisms (GMOs) of plants, animals or microorganisms, whose genetic material has been altered using genetic engineering techniques. However, GMOs are still controversial and raising some concerns over food safety in long terms .
The effective management of IPs mainly depends on chemical control methods so far, such as the application of pesticides, which is the quickest and most hard-hitting control method . Many pesticides are also involved in the enhancement of agricultural production through the expurgation of soil-borne pathogens. In paddy fields, nearly 15% of the total plant protection, chemicals are used for crop production . Among them, Abamectin (ABM) and Thiamethoxam (TXM) are the most impelling systematic pesticides widely used for rice, soybean, sunflower, cotton and potato seed treatments as well as in fields nowadays [6, 7].
ABM, the pesticide used to treat IPs, naturally generated as fermentation products by Streptomyces avermitilis, a soil actinomycete . ABM blocks nerve and muscle cells of the insects mostly by enhancing the effects of glutamate at the invertebrate-specific glutamate-gated chloride channel with minor impact on gamma-aminobutyric acid receptors [9,10,11]. This barricade causes an influx of chloride ions into the cells, leading to a hyper polarization and subsequent paralysis of invertebrate neuromuscular systems, while comparable doses are not toxic for mammals, as they do not possess glutamate-gated chloride channels .
TXM is a neonicotinoid that can be absorbed quickly by plants and transported to all of its tissues, including pollens where it acts to deter insect-feeding. This compound interferes with nicotinic acetylcholine receptors in the central nervous system of insects, and eventually paralyzes their muscular movements . TXM has been widely used because it controls a broad range of IPs while possessing relatively low mammalian toxicity [14, 15].
Although it is clear that pesticides can kill crop insects, it is still elusive whether they can affect plant growth and physiological performance . Generally, we could not see an obvious alteration of plant development after spraying a commercial pesticide, but this doesn’t mean that pesticide can’t influence the endogenous metabolic processes of crops, which may indirectly bring about human health issues. Thus, evaluating the effects of pesticides on crop physiology are crucial for IP control programs.
Rice (Oryza sativa L.) is a major staple cereal in the world, providing essential caloric requirements for more than half of the world’s population. To satisfy the gradually increasing food demands for a rapidly growing population, rice yields need to be increased up to 40% by 2030 . Meanwhile, many rice insects including brown plant hopper, leaf roller, and stem borer result in a major threat to rice production. To date, diverse insecticides have been used to suppress rice pests in an open field; among them, the application of ABM and TXM is the major solution for killing masticatory and sucking IPs. However, apart from crop safety, it is unknown whether plant physiology is compromised by the two pesticides, thereby triggering an interesting question to be addressed.
RNA sequencing (RNA-Seq) is a powerful tool to examine the continuously changing cellular transcriptome, thereby facilitates the ability to know potential physiological changes under distinct conditions. In this study, we conducted RNA-Seq analysis to determine rice dynamic performances after ABM and TXM spray through characterization of Differentially Expressed Genes (DEGs), Differentially Expressed Alternatively Spliced RNAs (DE AS), Differentially Expressed Long Non-Coding RNAs (DE lncRNAs) and Differentially Expressed Transposable elements (DE TEs). We found that a limited number of these coding and non-coding transcripts can be overlapped or exclusively changed along with the application of two different pesticides. These results provide valuable insights into the proper usage of pesticides against masticatory and sucking IPs in crops.
Identification of differentially expressed genes (DEGs) under Abamectin (ABM) treated rice
Pesticides can kill crop IPs, but their influence on different biological and physiological processes are still elusive. To investigate rice transcriptome in response to pesticides, we carried out RNA-Seq and measured FPKM values of genes under ABM treatment. A total of 470 DEGs were annotated in rice under ABM treatments. Correlation coefficients (R) of all the treatments were near to 1, showing a high repetition of the experiment in terms of data analysis, expression and sequence coverage (Fig. 1a, Additional file 1). To determine reliability in the transcriptome gene expressions (GE) profiles in ABM treatments, we randomly checked the expression patterns of six DEGs using RT-qPCR. Expression patterns of all the examined genes were similar to RNA-seq data, indicating the credibility of our transcriptome dataset for gene exploration (Additional file 2). Hence, it would be reliable to find out the influence of pesticide by our RNA-seq dataset. We compared DEGs with other expressed genes in relation to their percentages, and got the highest number of 192 DEGs (1.00%) under 1 day (1d) ABM treated plants, followed by 179 (0.91%) DEGs under 3 h (3 h) treatment, and 157 (0.83%) DEGs under 3 days (3d) treatment. These results indicated that DEGs were less in number compared to non-altered genes, and further implicated that the insecticide has a little grasp on GE level, mostly impacting 1d treated plants (Fig. 1b).
To further investigate the potential functions of DEGs, we identified their localization into different cellular components or biological processes under GO terms (Fig. 1c). Besides specifically expressed DEGs under three treatments of ABM, there were still some overlaps among DEGs per time point (Fig. 1d, Additional file 3), e.g., 3 h and 1d treatments shared 18 DEGs, six between 1d and 3d treatments, while 28 shared DEGs were recorded among 3d and 3 h. Apart from this, we also have three co-expressed DEGs shared by all treatments (Fig. 1d, Additional file 3). To further pursue dynamic changes in DEGs, we measured the FPKM values of genes under different time treatments of ABM. We observed two genes, Os01g69070 and Os06g45970, which are involved in Auxin response , were particularly induced, suggesting a potential alteration in auxin signaling by ABM treatments (Fig. 1e).
Identification of DEGs under Thiamethoxam (TXM) applied rice
Since we didn’t observe a severe alteration in GE after spraying ABM, we attempted to select another commercial pesticide to check whether it can lead to a significant change in rice transcriptomes. Due to the widespread use of TXM, it would be worthwhile to study its influence on the endogenous metabolic processes in plants.
TXM has a little more impact on GE level compared to ABM, as a total of 670 DEGs were detected in TXM treated rice. Reliability upon experiment was checked by correlation coefficients (R) which were near to 1 for all treatments (Fig. 2a). To further prevail the effectiveness of TXM, we compared DEGs with other expressed genes in terms of their percentages, and got highest number of DEGs 553 (2.94%) under 1d treatment of TXM, followed by expressions of 99 (0.52%) DEGs under 3d treatment, and 52 (0.27%) DEGs under 3 h treatments (Fig. 2b), adumbrating the fluctuating influence of this pesticide. DEGs were then annotated into functional categories using negative log10 (P-value), which illustrated their involvement into transport, localization or response to stimuli GO terms (Fig. 2c). Besides specifically expressed genes, there was a very small proportion overlap in DEGs per time point by TXM (Fig. 2d, Additional file 3). Furthermore, FPKM values of various DEGs with the abiotic stimulus, plastid, and transporter activity have dynamically changed in response to pesticide treatments (Fig. 2e), indicating that the application of TXM can induce some stress responses in rice.
Identification and characterization of the co-expressed DEGs by two pesticides
DEGs co-expressed among pesticides treatments are of prime importance due to their responses to both pesticides. After examining the individual effects of two pesticides, we scrutinized the co-expressed DEGs, and found 166 shared DEGs expressed under both insecticides treatments (Fig. 3a). To further study the localization and potency of the co-expressed DEGs, we carried out MapMan analysis in detail. These DEGs were mapped into hormone metabolism, RNA, stress, miscellaneous, protein, and signaling pathways with proportionate percentages of 6.62, 6.62, 6.02, 6.02, 4.81 and 3.61%, respectively (Fig. 3b). Notably, we found a significant upregulation of Os12g27220, which encodes Spermidine hydroxyl cinnamoyl transferase 1, an enzyme responsible for the biosynthesis of alkaloids, terpenoids, and phenolics (Fig. 3c) [19, 20], proclaiming more synthesis of spermidine may provide plants protection from diseases and pests by using agricultural chemicals.
Next, we examined DEGs specific to each drug by MapMan analysis. We found ABM-specific DEGs involved in several important processes, including response to stimuli, signaling, transport and protein (Fig. 3d, e). Interestingly, we noticed an induced expression level of some DEGs under 1d treatments of ABM compared to control or TXM, but decreased apparently at 3d treatment, indicating no longer effects of ABM on the expressions of these DEGs (Fig. 3f).
Specific DEGs under TXM treatments, by contrast, involved in different localization-related GO terms and cellular processes (Fig. 3g, h). Selected DEGs with unstable expressions indicated that TXM similar to ABM, has limited lasting effects on some GEs in rice (Fig. 3i, Additional file 3). Taken together, these data enlightened the limited roles of the two pesticides in GE regulation.
Identification of alternative splicing (AS) events in pesticides applied rice
In addition to be used for DEGs, RNA-seq dataset is also a good resource for AS analysis. Thus, we examined the AS changes affected by pesticides, and acquired approximately 3725 genes undergoing 5779 AS events (Additional file 4). Of them, 270 genes experienced 274 Differentially Expressed Alternative Splicing (DE AS) activity under both pesticides treatments (Fig. 4a and Additional file 5).
We classified total DE AS into four types; i.e., Exon skipping (ES), the most abundant (67.15%) of AS events, followed by Alternative 3′ splice site (A3SS) (20.44%), Intron retention (IR) (6.93%) and Alternative 5′ splice site (A5SS) (5.47%) (Fig. 4a and Additional file 5). We further investigated how AS contributed specifically to either ABM or TXM. For this purpose, we measured 178 DE AS events under ABM and 167 DE AS events under TXM, respectively, consisting of special and shared events among treatments (Fig. 4b, c and Additional file 5). Interestingly, we perceived the highest number of DE AS events under 1d treatments for both pesticides (special or shared), while this number reduced again under 3d treatment as in DEGs (Fig. 4d and Additional file 5), showing a dynamic alternation of AS triggered by these two agricultural chemicals. We predicted an A3SS AS event of Os03g60430, a AP2 domain protein-encoding gene, showing high AS activity under ABM compared with control samples under 1d treatment (Fig. 4e). Gene model presented a phenomenon that AS activity has happened in ABM treated samples, avoiding portion from the DEG on its 3′ site (Fig. 4e and f). These results evident that AS alterations could be occurred with pesticides.
Concurrently, we compared DEGs undergoing AS and transcriptional changes at the same time in response to studied pesticides, and got 11 DE AS events (Fig. 5a and b). As an illustration, the AS activity and gene expression of Os03g12620, which encodes Glycosyl hydrolases family 17, are oppositely regulated by pesticides (Fig. 5c).
Subsequently, we carried out the distribution of DEGs by DE AS into different pathways using MapMan analysis, and found them enriched in RNA and protein metabolic processes (Fig. 5d). An interesting example is that a gene model of the enzyme, Enoyl-CoA hydratase/isomerase family protein, declared the skipping of one exon under 1d ABM samples, accompanied by an increase of expression under ABM treatment compared to control, exclaiming AS involvement in the GE regulation after a drug spray (Fig. 5e and f). Further research will be interesting to explore the biological significance of DEGs by DE AS in plants under insecticidal environments.
Characterization of long non-coding RNAs (lncRNAs) in rice under pesticides treatments
LncRNA has been implicated playing a critical role in coding gene expressions. To predict lncRNAs in our transcriptomic dataset, we analyzed the assembled and filtered transcripts procuring approximately 3994 unique lncRNAs under two pesticides treatments, with 83 differentially expressed lncRNAs  among them (Fig. 6a). These differentially expressed lncRNAs (DELs) ranged in length between 270 and 6317 bp, and the most abundant length was 300–500 bp (Fig. 6a, Additional file 6). Furthermore, we distributed DELs individually to ABM and TXM treatments, and found the maximum number of DELs under TXM, consistent with the trend of DEGs, indicating a broader spectrum characteristics of TXM than ABM (Fig. 6b and c). We also examined co-expressed DELs among pesticides treatments and observed a higher number of DELs under 1d treatments (Fig. 6d), submitting a dynamic change of lncRNAs similar to DE AS in pesticides treated rice.
Previous studies have shown that lncRNAs could regulate the expressions of their neighboring protein-coding genes (PCGs) [22,23,24]. Therefore, we performed hierarchical clustering of the DELs, in which most of the DELs have positively regulated expressions of their neighboring genes (Fig. 6e). As an example, it is obvious that lncRNA TU37692 is positively regulating the expression of its neighboring Polysaccharide-K gene Os05g11260 (Fig. 6f).
A major theme involves, is the regulatory role of lncRNAs, which acts as a miRNA “sponge” to trigger the expression of PCGs [25,26,27]. To predict miRNA, lncRNA and coding genes interactions, we used Cytoscape (http://www.cytoscape.org/) and constructed the putative interactive network of miRNAs targeting their presumed lncRNAs and PCGs. Both, lncRNA TU9050A and Os08g37700 are predicted as miR1436 targets; we observed that they could be simultaneously induced by the two pesticides, suggesting that lncRNA TU9050A may block miR1436 activity to accelerate Os08g37700 transcription (Fig. 6g, h and Additional file 7). Likewise, Os01g55880, a target of miR2864.2, could be more induced through the sponge activity of lncRNA TU38959 under drugs treatments (Fig. 6h and Additional file 8). Nevertheless, we found that TU36112, the non-coding target of miR1858, and TU389592, the non-coding target of miR1846c-3p, could be dynamically changed, whereas the coding targets of miR1858 and miR1846c-3p were not profoundly altered by pesticide treatments. Thus, whether these two DELs can function as a miRNA target mimicry is unknown and should be further determined (Fig. 6i). Together, this predicted evidence suggest that lncRNAs might be involved in the fluctuating GE of DEGs through regulation of miRNA activities in response to pesticides.
Transposable elements (TEs) are involved in the alteration of gene expressions in the locale
TE insertion is another significant contributor of gene expressions in plants other than lncRNAs. In RNA-seq dataset, TE transcriptions can easily be examined as an asset. So we examined TEs in genes having transcriptional changes in response to pesticides treatments. Firstly, we analyzed DE TEs individually for two pesticides and got 193 and 387 DE TEs under ABM and TXM treatments, respectively (Fig. 7a, b and c). Intriguingly, DE TEs were higher in TXM treatments as that of DE AS, DEGs as well as DELs, further illustrating TXM has stronger effects than ABM on rice. Ninety-two co-expressed DE TEs were ascertained responsive to both pesticides (Fig. 7c, Additional file 9). Furthermore, we classified DE TEs on the basis of their functions and observed that the maximum number of DE TEs are in the class Miniature Inverted-repeat Transposable Elements 175 (35.9%), followed by retrotransposon 164 (33.6%) and transposons 149 (30.5%) (Fig. 7d).
The density of TEs to the nearest genes can result in impacts on transposition conversely affecting transcription. Therefore, we measured the density of DE TEs and TEs to the nearest genes indicated by a peak (P-value < 2.2e-16) and found that DE TEs were closer to genes compared to random TEs (Fig. 7e). These results indicated that DE TEs have a higher ability to regulate the nearby genes.
Besides coding regions, GE levels can also be impacted by TE insertions into non-coding regions. We found that more than 42% of DE TEs were overlapped with lncRNAs in our RNA-seq datasets, confirming the phenomenon that GE regulation by TEs can be the result of its insertion into either coding or non-coding DNA zones (Fig. 7f). Furthermore, hierarchical interaction showed that most of the DE TEs were involved in positively regulating their mediated GE levels, for instance, TE-48405 on rice chromosome 12 near to the TSS of the Os12g38180, a gene encoding Heat shock cognate 70–1, positively regulated its expression pattern in the vicinity by ABM and TXM treatments (Fig. 7g, h and Additional file 9). The above results suggested that TEs are the crucial means to alter gene transcriptions, playing critical roles in stress responses and defense mechanisms.
As a plant development hormone, auxins play critical roles in plant growth and reproduction. In the present study, we found that ABM regulates many genes with potential functions in the auxin response, including one auxin efflux carrier gene Os01g69070, three auxin-repressed genes Os03g22270, Os11g44810, and Os08g35190, and two auxin response SAUR genes Os06g45970, Os02g07110, suggesting that ABM may influence the auxin signaling in rice (Fig. 1e). Under stressful conditions, plants would produce a high level of anthocyanins to increase yield and antioxidant capacity . We found that the expression of Os04g54550 (AT5MAT) encoding an enzyme that prevents the degradation of anthocyanins, can be induced considerably by TXM (Fig. 2e). Hence, upregulation of AT5MAT may be an adaptive strategy of plants to cope with the toxic effects of insecticide . Heat Shock Proteins (HSPs) are involved in the regulation of specific substrate proteins in stressful conditions especially under high temperature . Under TXM application, some HSPs were also inhibited (Fig. 2e), showing the possibility that disruption of thermal morphogenesis by the farm chemical is an alternative approach to trade off the physiological and defense responses in plants.
Thaumatins are pathogenesis-related (PR) proteins, which are induced by various agents ranging from ethylene to pathogens in plants . Therefore, the identification and functional characterization of thaumatins were pretty interesting. We found that Os03g46070 and Os12g43380, two genes involved in thaumatins biogenesis were increased under ABM treatments (Fig. 3f), indicating thaumatins may play coordinated roles in defense against masticatory insects. Nonetheless, salicylic acid (SA) and jasmonic acid (JA), two phytohormones play roles in resistance to pathogens by inducing the production of pathogenesis-related proteins and are further involved in the systematic acquired resistance of plants where a pathogenic attack on one part of the plant induces resistance in other parts [32, 33]. Here, two genes were alternatively regulated for SA production inducing the expressions of Os05g01140, while decrease in the expressions of Os11g15040 were illustrated under two pesticides treatments. Furthermore, Os09g02710, Os08g39850, Os12g37260 and Os03g28940, four genes involved in JA synthesis were increased during pesticides application (Additional file 10). These results suggested that application of studied pesticides reprogramed different metabolic pathways by modulating the expressions of critical DEGs in rice.
Os11g32880, a DEAD-box ATP-dependent RNA helicase (DBH) encoding gene, was subjected to AS and transcriptionally regulated by TXM treatments (Fig. 5b). It is noteworthy that there are a lot of redundant helicases in rice, which might be the result of either tandem or segmental duplications , thus whether AS of Os11g32880 can regulate other helicases by pesticides remains unknown, since AS events sometimes exert dominate-negative effects on gene functions . As a defense signal to IPs, terpene synthase, which is involved in secondary metabolism and biosynthesis of phytoalexins, was negatively regulated by high ES under both ABM and TXM treatments , suggesting that pesticide can regulate these defensive genes at post-transcriptional level to enhance plant biotic resistance. Taken together, the subsequent induction in expression of many stress-responsive genes indicated that AS might be involved to persuade defense mechanism of rice in response to pesticides.
In a meanwhile, this study was motivated by our interest in the spatio-insecticidal incentives of lncRNAs, and exploring its potential biological roles in GE modulation in rice treated with pesticides. Interestingly, we found an elevated number of DELs at 1d treatment under both pesticides, which evidenced that the increased number of DEGs might be due to the involvement of adjacent lncRNAs (Fig. 6b and c). Add-on to this, most of the lncRNA-mediated genes belong to either metabolic, stress-responsive pathways or involved in defense mechanisms (Fig. 6e). Through the potential to create new transcripts, premature termination, or novel promoters, TEs can influence the expression patterns of plants genes [37, 38]. Here, we have evaluated a variety of TEs that can interact with neighboring genes to shape transcripts and regulate their expression levels. For instance, C-glycosylflavone coding gene, Os06g01250, is critical for defense against plant-eating insects [39,40,41], was positively regulated by TE-234312 in rice under both chemicals, determining its roles in defense response probably via TE mediated regulations (Fig. 7g). NB-ARC domain is found in bacteria and eukaryotes and shared by plant resistance (R) proteins proposed to be involved in pathogen recognition . Interestingly, we found Os08g09110; an NB-ARC R gene was negatively regulated by TE-323079 under both pesticides treatments, indicating agricultural drugs may affect R gene expressions indirectly by regulating TEs.
To conclude, this study provides the first systematic analysis of gene expressions responsive to two commercial pesticides in rice, which can be further applied for other crops and alternative agricultural chemicals. Although we detected several DE AS, DELs and DE TEs alterations accompanied with host GEs, our study overall depicted that recommended doses of ABM and TXM pesticides are effective in controlling IPs and likely to have no apparently harmful effects on plant growth and physiological performance in rice.
Rice (Oryza sativa), cv. Xidao No.1, spp. japonica WT was used in the study. Seeds were soaked in deionized water overnight at 28οC and then transferred to a net floating on a 0.5 mM CaCl2 solution. After 7 days, seedlings were transplanted into soil in a greenhouse (approximately 28/25οC day/night, 75 ± 85% relative humidity) and grown for 6 weeks. Pesticides ABM and TXM were applied at the jointing stage (6 weeks) using the recommended rate; i.e.,1.81 gh− 1 of ABM and 3 g ai 100 l− 1, (300 l h− 1) of TXM . Pesticides were applied at 8 o’clock in the morning and samples (ABM, TXM and control) were collected at three intervals; i.e., 3 h, 1 day and 3 days after pesticides spray. All nine samples (each one containing three biological replicates) were immediately frozen in liquid nitrogen and stored at − 80 οC for RNA extraction.
RNA extraction and Illumina sequencing
Three biological replicates were used for all RNA-Seq experiments sampled from control and pesticide-treated plants. Total RNA was extracted from each sample using TRIzol Reagent (Invitrogen)/RNeasy Mini Kit (Qiagen). Total RNA was quantified and qualified by Agilent 2100 Bio-analyzer (Agilent Technologies, Palo Alto, CA, USA), Nano Drop (Thermo Fisher Scientific Inc.) and 1% agarose gel. RNA with RIN value above 7 was used for further sequencing library construction. RNA-sequencing was performed on an Illumina HiSeq2500 platform (PE 150 bp) at Novogene Company (Beijing, China). The detailed information for RNA-seq data is listed in supporting information Additional file 1. All the RNA-seq data have been deposited in NCBI database with accession number PRJNA532802.
Analysis of RNA-seq data
Rice genomic sequence (O. sativa_323_v7.0) and the corresponding annotation was retrieved from the Phytozome database (Version 12.0) . TEs were extracted from the rice repeat annotation file (Osativa_323_v7.0.repeatmasked_assembly_v7.0.gff3). Only TEs derived from the intergenic region remained for further analysis. After filtering adaptor and low quality reads with FASTX-Toolkit (version 0.0.14), all clean data for each sample were separately aligned to the reference genome using Hisat2 with the following parameters: “--end-to-end, --known-splice site-infile” . For each sample, unique mapped reads were assembled into putative transcripts based on a reference-guided assembly strategy using StringTie (version 1.3.3b) . The meta-assembly tool TACO (version 0.7.3) was used to merge putative transcripts from each sample into a unified set of transcripts , which was then compared to the reference gene GTF file using gffcompare (version 0.10.1) .
Based on transcript classification codes, “u” (unknown) intergenic transcripts were regarded as novel gene loci and used to lncRNA identification. Five steps were adopted to identify bona fide lncRNAs as previously described : (1) transcripts should be with length ≥ 200 bp and detected in more than 3 samples; (2) transcripts derived from rRNA and tRNA were removed (cutoff E-value0.001); (3) transcripts encoding proteins and protein-coding domains were removed by searched against the Swiss-Prot and Pfam databases (cutoff E-value 0.001); (4) OrfPredictor was applied to predict ORFs and transcripts that encode more than 100 amino acid was removed . (5) transcripts were removed that did not pass the protein-coding-score test using Coding-Non-Coding Index (CNCI)  and Coding Potential Calculator (CPC) software .
Resulting lncRNA transcripts and known transcripts, intergenic TEs were then merged into non-redundant transcripts, which were further quantified by StringTie for each sample . Differential expression analysis for each sequenced library was performed using ballgown . The corrected P value of 0.05 and abs |log2 (Fold change)| of 1 were set as the threshold for significant differential expression. Singular Enrichment Analysis from AgriGO was performed to identify significantly enriched GO terms in the gene list out of the background of the reference gene list . GO terms pathways with false discovery rate (q-value) < 0.05 were considered as significantly altered.
All putative AS events were extracted from the rice transcript GTF file using rMATS . Expressed AS events were identified by filtered below three samples with total IJC + SJC < 30. Differentially AS events between control and each pesticide-treated condition were identified using rMATS with a stringent threshold: FDR ≤ 0.05, Delta PSI (ΔPSI) ≥10%.
Prediction of miRNA, lncRNA and coding genes network
Pearson correlation was employed to explore the expression relationship between lncRNAs and their neighboring genes (≤ 10 Kb). Mature miRNAs for rice were retrieved from the miRBase database (Release 21.0). To identify lncRNAs and genes as putative miRNA targets and mimics, the TAPIR tool was used with the default settings . The relationship between miRNAs, lncRNAs, and genes were used to construct the interaction networks with Cytoscape software (version 3.5.1) .
Validation of DEGs by RT-qPCR
To validate RNA-seq, we selected six genes based on their expressions from DEGs and tested their expressions using qRT-PCR. RNA used earlier for RNA-seq samples were reversely transcribed using reverse transcriptase cDNA synthesis kit (GoScript™ Reverse Transcription System, Promega, Beijing, China). qRT-PCR analyses of genes were performed as previously described . Real-time qPCR was done in triplicates on Step-One Plus real-time PCR system (ABI) with the Power Up SYBR Master Mix (ABI) with three biological replicates. The following cycling conditions were used for Real-time qPCR: 2 min at 95 °C, 40 cycles of 10 s at 95 °C and 40 s at 65 °C, and a final step for melting curve determination (15 s at 95 °C, 1 min at 60 °C and 15 s at 95 °C). ACTIN was used as an internal control. Gene expression was calculated using the 2-ΔΔCt method. Primers used for gene expression analysis are listed in Additional file 11.
Availability of data and materials
All the raw reads produced in this study have been deposited in NCBI database with accession number PRJNA532802.
Alternative 3′ Splice Site
Alternative 5′ Splice Site
Coding Potential Calculator
- DE AS:
Differetially Expressed Alternatively Spliced
- DE TEs:
Differentially Expressed Transposable Elements
Differentially Expressed Genes
Differentially Expressed Long Non-coding RNAs
Fragments Per Kilobase of exon per Million reads
Long Non-conding Ribonucleic Acid
Mutually Exclusive Exon
Protein Coding Genes
Bardner R, Fletcher K. Insect infestations and their effects on the growth and yield of field crops: a review. Bull Entomol Res. 1974;64(1):141–60.
Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, Naylor RL. Increase in crop losses to insect pests in a warming climate. Science. 2018;361(6405):916–9.
Lombardo L. Genetic use restriction technologies: a review. Plant Biotechnol J. 2014;12(8):995–1005.
Mamta B, Rajam M. RNAi technology: a new platform for crop pest control. Physiol Mol Biol Plants. 2017;23(3):487–501.
Agnihotri N. Pesticide consumption in agriculture in India-an update. Pestic Res J. 2000;12(1):150–5.
Alford A, Krupke CH. Translocation of the neonicotinoid seed treatment clothianidin in maize. PLoS One. 2017;12(3):e0173836.
Cheng X, Wang Y, Li W, Li Q, Luo P, Ye Q. Nonstereoselective foliar absorption and translocation of cycloxaprid, a novel chiral neonicotinoid, in Chinese cabbage. Environ Pollut. 2019;252:1593–8.
Ikeda H, Nonomiya T, Usami M, Ohta T, Ōmura S. Organization of the biosynthetic gene cluster for the polyketide anthelmintic macrolide avermectin in Streptomyces avermitilis. Proc Natl Acad Sci. 1999;96(17):9509–14.
Bloomquist JR. Chloride channels as tools for developing selective insecticides. Arch Insect Biochem Physiol. 2003;54(4):145–56.
Cully DF, Vassilatis DK, Liu KK, Paress PS, Van der Ploeg LH, Schaeffer JM, Arena JP. Cloning of an avermectin-sensitive glutamate-gated chloride channel from Caenorhabditis elegans. Nature. 1994;371(6499):707.
Bloomquist JR. Ion channels as targets for insecticides. Annu Rev Entomol. 1996;41(1):163–90.
de Faria DBG, Montalvão MF, de Souza JM, de Oliveira MB, Malafaia G, de Lima Rodrigues AS. Analysis of various effects of abamectin on erythrocyte morphology in Japanese quails (Coturnix japonica). Environ Sci Pollut Res. 2018;25(3):2450–6.
Mörtl M, Kereki O, Darvas B, Klátyik S, Vehovszky Á, Győri J, Székács A: Study on soil mobility of two neonicotinoid insecticides. J Chem. 2016. p. 4546584.
Maienfisch P, Angst M, Brandl F, Fischer W, Hofer D, Kayser H, Kobel W, Rindlisbacher A, Senn R, Steinemann A. Chemistry and biology of thiamethoxam: a second generation neonicotinoid. Pest Manag Sci. 2001;57(10):906–13.
Reisig DD, Herbert DA, Malone S. Impact of neonicotinoid seed treatments on thrips (Thysanoptera: Thripidae) and soybean yield in Virginia and North Carolina. J Econ Entomol. 2012;105(3):884–9.
Afifi M, Lee E, Lukens L, Swanton C. Thiamethoxam as a seed treatment alters the physiological response of maize (Zea mays) seedlings to neighbouring weeds. Pest Manag Sci. 2015;71(4):505–14.
Khush GS. What it will take to feed 5.0 billion rice consumers in 2030. Plant Mol Biol. 2005;59(1):1–6.
Wang J-R, Hu H, Wang G-H, Li J, Chen J-Y, Wu P. Expression of PIN genes in rice (Oryza sativa L.): tissue specificity and regulation by hormones. Mol Plant. 2009;2(4):823–31.
Campos L, Lisón P, López-Gresa MP, Rodrigo I, Zacarés L, Conejero V, Bellés JM. Transgenic tomato plants overexpressing tyramine N-hydroxycinnamoyltransferase exhibit elevated hydroxycinnamic acid amide levels and enhanced resistance to Pseudomonas syringae. Mol Plant-Microbe Interact. 2014;27(10):1159–69.
Walters D, Meurer-Grimes B, Rovira I. Antifungal activity of three spermidine conjugates. FEMS Microbiol Lett. 2001;201(2):255–8.
Si-Ammour A, Windels D, Arn-Bouldoires E, Kutter C, Ailhas J, Meins F, Vazquez F. miR393 and secondary siRNAs regulate expression of the TIR1/AFB2 auxin receptor clade and auxin-related development of Arabidopsis leaves. Plant Physiol. 2011;157(2):683–91.
Wang T-Z, Liu M, Zhao M-G, Chen R, Zhang W-H. Identification and characterization of long non-coding RNAs involved in osmotic and salt stress in Medicago truncatula using genome-wide high-throughput sequencing. BMC Plant Biol. 2015;15(1):131.
Liu J, Wang H, Chua NH. Long noncoding RNA transcriptome of plants. Plant Biotechnol J. 2015;13(3):319–28.
Deng P, Liu S, Nie X, Weining S, Wu L. Conservation analysis of long non-coding RNAs in plants. Sci China Life Sci. 2018;61(2):190–8.
Dykes IM, Emanueli C. Transcriptional and post-transcriptional gene regulation by long non-coding RNA. Genomics Proteomics Bioinformatics. 2017;15(3):177–86.
Wu HJ, Wang ZM, Wang M, Wang XJ. Widespread long noncoding RNAs as endogenous target mimics for MicroRNAs in plants. Plant Physiol. 2013;161(4):1875–84.
Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, Leyva A, Weigel D, Garcia JA, Paz-Ares J. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39(8):1033–7.
Kovinich N, Kayanja G, Chanoca A, Riedl K, Otegui MS, Grotewold E. Not all anthocyanins are born equal: distinct patterns induced by stress in Arabidopsis. Planta. 2014;240(5):931–40.
Wang H, Fan W, Li H, Yang J, Huang J, Zhang P. Functional characterization of dihydroflavonol-4-reductase in anthocyanin biosynthesis of purple sweet potato underlies the direct evidence of anthocyanins function against abiotic stresses. PLoS One. 2013;8(11):e78484.
Bernfur K, Rutsdottir G, Emanuelsson C. The chloroplast-localized small heat shock protein Hsp21 associates with the thylakoid membranes in heat-stressed plants. Protein Sci. 2017;26(9):1773–84.
Ruiz-Medrano R, Jimenez-Moraila B, Herrera-Estrella L, Rivera-Bustamante RF. Nucleotide sequence of an osmotin-like cDNA induced in tomato during viroid infection. Plant Mol Biol. 1992;20(6):1199–202.
Hayat S, Ali B, Ahmad A. Salicylic acid: biosynthesis, metabolism and physiological role in plants. In: Hayat S, Ahmad A (eds) Salicylic acid: A plant hormone. Drodrecht: Springer; 2007. p. 1–14.
Lyons R, Manners JM, Kazan K. Jasmonate biosynthesis and signaling in monocots: a comparative overview. Plant Cell Rep. 2013;32(6):815–27.
Yu J, Wang J, Lin W, Li S, Li H, Zhou J, Ni P, Dong W, Hu S, Zeng C. The genomes of Oryza sativa: a history of duplications. PLoS Biol. 2005;3(2):e38.
Qin Z, Wu J, Geng S, Feng N, Chen F, Kong X, Song G, Chen K, Li A, Mao L, et al. Regulation of FT splicing by an endogenous cue in temperate grasses. Nat Commun. 2017;8:14320.
Li R, Tee C-S, Jiang Y-L, Jiang X-Y, Venkatesh PN, Sarojam R, Ye J. A terpenoid phytoalexin plays a role in basal defense of Nicotiana benthamiana against potato virus X. Sci Rep. 2015;5:9682.
Elbarbary RA, Lucas BA, Maquat LE. Retrotransposons as regulators of gene expression. Science. 2016;351(6274):aac7247.
Lisch D. How important are transposons for plant evolution? Nat Rev Genet. 2013;14(1):49.
Ma Y, Li H, Lin B, Wang G, Qin M. C-glycosylflavones from the leaves of Iris tectorum maxim. Acta Pharm Sin B. 2012;2(6):598–601.
Zhou G, Wu H, Wang T, Guo R, Xu J, Zhang Q, Tang L, Wang Z. C-glycosylflavone with rotational isomers from Vaccaria hispanica (miller) Rauschert seeds. Phytochem Lett. 2017;19:241–7.
Du Y, Chu H, Chu IK, Lo C. CYP93G2 is a flavanone 2-hydroxylase required for C-glycosylflavone biosynthesis in rice. Plant Physiol. 2010;154(1):324–33.
Uehling J, Deveau A, Paoletti M. Do fungi have an innate immune response? An NLR-based comparison to plant and animal immune systems. PLoS Pathog. 2017;13(10):e1006578.
Nordberg H, Cantor M, Dusheyko S, Hua S, Poliakov A, Shabalov I, Smirnova T, Grigoriev IV, Dubchak I. The genome portal of the Department of Energy Joint Genome Institute: 2014 updates. Nucleic Acids Res. 2013;42(D1):D26–31.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650.
Niknafs YS, Pandian B, Iyer HK, Chinnaiyan AM, Iyer MK. TACO produces robust multisample transcriptome assemblies from RNA-seq. Nat Methods. 2016;14(1):68.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7(3):562.
Wang M, Yuan D, Tu L, Gao W, He Y, Hu H, Wang P, Liu N, Lindsey K, Zhang X. Long noncoding RNAs and their proposed functions in fibre development of cotton (Gossypium spp.). New Phytol. 2015;207(4):1181–97.
Min XJ, Butler G, Storms R, Tsang A. OrfPredictor: predicting protein-coding regions in EST-derived sequences. Nucleic Acids Res. 2005;33(suppl_2):W677–80.
Mukherjee S, Manna S, Mukherjee P, Panda CK. Differential alterations in metabolic pattern of the spliceosomal uridylic acid-rich small nuclear RNAs (UsnRNAs) during malignant transformation of 20-methylcholanthrene-induced mouse CNCI-PM-20 embryonic fibroblasts. Mol Carcinog. 2009;48(9):773–8.
Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(suppl_2):W345–9.
Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(suppl_2):W64–70.
Shen S, Park JW, Lu Z-X, Lin L, Henry MD, Wu YN, Zhou Q, Xing Y. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci. 2014;111(51):E5593–601.
Bonnet E, He Y, Billiau K, Van de Peer Y. TAPIR, a web server for the prediction of plant microRNA targets, including target mimics. Bioinformatics. 2010;26(12):1566–8.
Saito R, Smoot ME, Ono K, Ruscheinski J, Wang P-L, Lotia S, Pico AR, Bader GD, Ideker T. A travel guide to Cytoscape plugins. Nat Methods. 2012;9(11):1069.
Qin Z, Bai Y, Muhammad S, Wu X, Deng P, Wu J, An H, Wu L. Divergent roles of FT-like 9 in flowering transition under different day lengths in Brachypodium distachyon. Nat Commun. 2019;10(1):812.
We thank Dr. Qingyao Shu for the gift of Xidao No.1 rice seeds.
In this study, the design of the study, RNA-sequencing and writing the manuscript were mainly supported by Zhejiang Provincial Natural Science Foundation of China (LR16C060001). Sample collection, data analysis and experimental validation were supported in part by the National Natural Science Foundation of China (91640109, 31871588), China Postdoctoral Science Foundation (2018 M630679), and Hundred-Talent Program of Zhejiang University.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Mapping statistics of mock and pesticides applied three treatments. (A-Abamectin; T-Thiamethoxam)
Confirmation of the expression patterns of DEGs using real-time quantitative polymerase chain reaction (RT-qPCR)
DEGs identified by RNA-seq analysis under two pesticides treatments expressed in FPKM
Pie chart represents all expressed genes and total AS activities observed in the study. All AS activities were divided into five types; i.e., Exon skipping (ES), Alternative 3′ splice site (A3SS), Alternative 5′ splice site (A5SS), Mutually exclusive exon (MXE) and Intron retention (IR). Percentage of each type is also acknowledged
DEGs identified executing AS events under two pesticides treatments
DELs identified under two pesticides treatments expressed in FPKM
Predicted interaction network of miRNAs, lncRNAs, and PCGs. Circles show PCGs, triangles represent miRNAs, and hexagonal structures indicate lncRNAs
Differentially expressed genes targeted by miRNAs under pesticides treatments
DE TEs identified in rice under two pesticides treatments
DEGs related to salicylic acid and jasmonic acid synthesis for rice under two pesticides treatments
List of Primers used for validation of six auxin responsive genes
About this article
Cite this article
Muhammad, S., Tan, J., Deng, P. et al. Pesticide application has little influence on coding and non-coding gene expressions in rice. BMC Genomics 20, 1009 (2019). https://doi.org/10.1186/s12864-019-6381-y
- Long non-coding RNAs
- Alternative splicing
- Abiotic stress