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 [18], 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 [21] 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.