Pesticide application has little influence on coding and non-coding gene expressions in rice

Abstract

Background: 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.

Results: 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.

Conclusion: 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.

Keywords: Abiotic stress; Alternative splicing; Long non-coding RNAs; Pesticide; Rice.

Fig. 1

Expression pattern and functional analysis of differentially expressed genes (DEGs) in rice inoculated with Abamectin (ABM). a Bar graphs depict correlation co-efficients (R) of ABM under three treatments, i.e., 3 h, 1d, and 3d. The y-axis represents correlation co-efficient of treatments, and x-axis shows pesticide treatments. b Proportionate percentages of DEGs to other expressed genes, red color in the bar graph shows the proportion of DEGs to other expressed genes illustrated in blue color. c Overview of Gene Ontology analysis of all DEGs under ABM application. The x-axis represents the negative log of the P-value, and y-axis shows GO terms. d Venn diagram describing total, unique and overlaps among DEGs after three treatments of ABM, the number of shared DEGs are specified in circles. e Expressions of selected DEGs based on high throughput sequencing, under control and ABM, treated plants. The y-axis is the FPKM (Fragments Per Kilobase of exon per Million reads) values for each gene and x-axis represents treatments of ABM. First two genes are the typical examples of induced genes under ABM compared with control, while others are examples for low expressed genes under ABM treatments

Fig. 2

Expression pattern and functional analysis of DEGs in rice inoculated with Thiamethoxam (TXM). a Bar graphs show correlation co-efficients (R) of three TXM treatments, i.e., 3 h, 1d, and 3d. The y-axis represents correlation co-efficient of treatments, and x-axis shows pesticide treatments. b Bar graphs represent proportionate percentages of DEGs to other expressed genes. The red color in the graph shows the proportion of DEGs to other expressed genes presented in blue color. c Overview of GO analysis of the putative DEGs under TXM application. The x-axis represents the negative logarithm of the P-value, and y-axis shows GO terms. d Venn diagram is describing total, unique and overlaps among DEGs after treatments with TXM. e Expressions of selected DEGs based on high throughput sequencing, under control and TXM treated plants. Expression levels in FPKM of the genes are given on y-axis along with their treatments on x-axis. The first three genes are typical examples of TXM which accumulate more under TXM treatments compared to control, while others are examples of low expressed genes under TXM treatments

Fig. 3

Expression profiles and functional distribution of co-expressed DEGs under ABM and TXM treatments. a Comparison of shared and unique DEGs under ABM and TXM treatments in rice. b MapMan pathway analysis for all co-expressed DEGs identified between ABM and TXM. The y-axis shows the distribution of genes into different pathways, while x-axis represents a number of genes assumed for each category. c The expression level of representative shared DEGs under control, ABM or TXM treatments. FPKM values are specified on y-axis, while x-axis represents treatment time. d Enriched GO terms of DEGs annotated in biological processes specific to ABM treatment. e Enriched MapMan pathways analyses for all unique DEGs expressed under ABM insecticide application. f DEGs specifically responsive to ABM treatments at different time intervals. g GO enrichment of TXM special DEGs. h Enriched MapMan pathways analyses for all unique DEGs of TXM. i Expression profiles of TXM representative DEGs under control, ABM or TXM treatments at three intervals. The expression level of genes is in FPKM, specified on y-axis, while x-axis represents treatment time

Fig. 4

Statistics of differentially expressed alternatively spliced (DE AS) genes under ABM and TXM treatments. a Pie chart represents all expressed DEGs (270) with Alternative Splicing (AS) events (274). AS genes along with their percentages are divided into four sub-categories; Exon skipping (ES), Alternative 3’splice site (A3SS), Alternative 5’splice site (A5SS), and Intron retention (IR). b Venn diagram represents shared and unique DEGs and their approximate DE AS under ABM treatments. A total number of DE AS (inside) and genes expressed (outside) at each treatment of ABM are specified along with treatment information. c Venn diagram represents shared or unique DEGs and their approximate DE AS under TXM treatments. d Graphical distribution of DE AS in response to ABM or TXM treatments. Total unique and shared DE AS concerning time are provided in squares or alongside arrows, respectively. e An example of A3SS of Os03g60430, AP2 domain-containing protein at relatively 1d treatment under control or ABM. Graphical representation of the gene showing AS activity at 3 site under ABM treated samples. f AS score (lncLevel) is predicted as an example under control, ABM and TXM treatments at three intervals. The y-axis shows the AS score, the highest is 1, while x-axis demonstrates the three treatments. Graph shows high AS activity of Os03g60430 under control and TXM at 1d treatment

Fig. 5

Relative expression and functional distribution of shared DEGs with DE AS under ABM and TXM treatments. a Venn diagram represents shared and unique DEGs and DE AS under ABM and TXM treatments. b Heat map represents the expression level of selected genes along with their functions and AS activity under insecticides treatments. Transcript levels following insecticides treatments are depicted using FPKM values on a color scale. The spots highlighted in Pink-magenta indicated the DEGs exhibit a significant expression level compared with control after treatments. c Expression level of the representative DE AS gene Os03g12620 under control, ABM and TXM treatments. The left side figure represents the AS score, while right side shows the expression level of AS-mediated gene depicting negative relationship. d Enriched MapMan pathways for DE AS events expressed under control, ABM and TXM treatments. The y-axis represents the distribution of genes into different cellular components, and x-axis shows gene numbers indicated in front of each bar. e An example of Exon skipping (ES) of Os06g39344 gene along with its AS score (lnc level) at 1d treatment under control or ABM treatments. Predicted graphical representation of the AS activity can be observed in the form of one exon skipping from the ABM treated samples. f Expression profile of representative DE AS gene under control, ABM and TXM treatments at three different intervals. The left side graph represents the AS score of DE AS, while right side shows the expression level of the gene with negative correlation reducing the expression level of the gene under pesticides treatments compared to control samples

Fig. 6

Expression profiles of differentially expressed long non-coding RNAs in rice exposed to two pesticides. a Line graph represents a total number of predicted expressed lncRNAs and protein-coding genes (PCG). Predicted length (aa) is shown on x-axis with scale, and cumulative frequency is revealed on y-axis. b Venn diagram shows shared and unique DELs perceived under ABM treatments. c Venn diagram represents DELs observed at each treatment of TXM. d Distribution of total, unique and shared DELs under three treatments responding ABM or TXM. Number of unique and shared lncRNAs are specified in squares or alongside arrows, respectively. e Heat map represents the expression level of lncRNAs and their mediated genes in response to the studied pesticides. Color scale indicates FPKM change (blue, low expression level and red, high expression level). Correlation specificity score is presented on the right side of the heat map for lncRNAs and its neighboring genes. Values close to 1 means high correlation (R) of DELs and genes in the vicinity. f The expression level of a selected lncRNA and its adjoining gene under ABM and TXM pesticides. Expression levels are assumed on y-axis and treatment time is denoted by x-axis. As an example, lncRNA TU37692 positive correlates Os05g11260. Graphical representation of the gene and lncRNA is shown above graph representing the position of gene and lncRNA. g Predicted interaction network of miRNAs, lncRNAs, and PCGs. Circles show PCGs, triangles represent miRNAs and hexagonal structures indicate lncRNAs. Osa-miR1436 is specified as an example, targeting lncRNA TU9050 and a gene Os08g37700 highlighted in the red color. h Expression levels of two lncRNAs-mediated genes targeted by miRNAs are specified. The y-axis represents expression levels of lncRNAs and PCG and x-axis shows treatments of two pesticides. Osa-miR1436 and Osa-miR2864.2 target lncRNAs and PCGs positively regulating their expression levels under drugs treatments. i Predicted base pair interaction between two miRNAs and their targeted lncRNAs. Line graphs show high expression level of lncRNAs under pesticide treatments compared with control

Fig. 7

Classification and expression profile of Transposable Elements (TEs) under ABM or TXM treatments. a Venn diagram represents shared, and unique DE TEs observed under three treatments of ABM. Total, unique (outside) and shared (inside) DE TEs are specified. b Shared and unique DE TEs perceived under different treatments of TXM. c Venn diagram depicts special, and co-expressed DE TEs expressed under ABM and TXM applications. d Bar graphs represent distribution of total DE TEs into three classes, i.e., Transposons 149, Retrotransposons 164 and Miniature Inverted-repeat TEs (MITE) 175 on the basis of their function. The y-axis represents number of expected DE TEs along with their percentages at the top of each bar and their classification is shown on x-axis. e Graph presents density of significant DE TEs and random TEs. The y-axis represents density of DE TEs and TEs (maximum value is 1) while x-axis shows distance to nearest genes (scale given in bp). P-value of significance is shown on the top of y-axis. f Pie chart represents proportions of DE TEs under pesticides treatments and interaction of DE TEs with DELs or lncRNAs. g Heat map shows FPKM values of TEs and adjoining genes in response to the studied pesticides. Color scale indicates expression level of TEs and their neighboring genes. Line graph on the right side of heat map represent correlation co-efficient (R) of TEs and DEGs. Most of the TEs positively correlates their neighboring genes. h Example of TE affecting expression level of its neighboring gene. Map shows graphical representation of TE 48405 and HS cognate protein 70–1 gene along with UTR and CDS regions indicating positive correlation of the TE expression and GE