Our main aim was to establish a resource to study transcriptomic patterns in A. myosuroides plants resistant or sensitive to leaf-applied herbicides in the absence of herbicide and at the early stages of response to herbicide, using experimental conditions as similar as possible to realistic field conditions. Herbicide damage to plants starts occurring 3 to 8 h after herbicide application (reviewed in [10]). To be efficient, NTSR must be implemented before herbicide damage is irreversible, and must be upheld long enough to allow resistant plants to recover [10]. Accordingly, RNA-Seq data was obtained from a time course experiment ranging from UT until 73HAT. A commercial iodosulfuron + mesosulfuron formulation that is applied in the field was used together with its recommended adjuvant in our experiments, because NTSR is selected for in the field by recurrent applications of not only herbicide molecules, but also associated formulations and adjuvants. Last, resistant and sensitive F2 plants derived from a pairing between a sensitive plant and a resistant plant from a population where NTSR had evolved under herbicide selective pressure were used as starting plant material because using plants with a similar genetic background was expected to facilitate identification of transcriptomic differences related to NTSR [10].
As there is no genomic resource for A. myosuroides, the Illumina technology was selected for sequencing because it is the technology of choice for de novo transcriptome deep sequencing and assembly without a reference genome [25]. After discarding assembly waste, ALOMYbase contained 65,558 contigs potentially from expressed genes that encoded 32,138 putative peptides (Table 1). This is to be compared to the 95 nucleotide sequences from A. myosuroides that had been deposited in GenBank/EMBL to the 5th of June 2015. Our work therefore tremendously increased the sequence data available for A. myosuroides.
The contigs in ALOMYbase were assembled for RNA-Seq data obtained from the aerial part of six F2 plants at the vegetative growth stage. Thus, ALOMYbase was not expected to be a comprehensive A. myosuroides transcriptomic resource. Accordingly, the completeness of ALOMYbase as estimated by CEGMA analysis was 61.3 %. Considering our starting material and CEGMA estimates of completeness for the transcript sets associated to the five grass genomes (79.8 to 99.2 %), A. myosuroides transcriptome representation in ALOMYbase was deemed satisfactory for our purpose.
As A. myosuroides genetic variability is high [40], the two parental plants used to generate the F2 population were most likely genetically different. Thus, nucleotide variation was certainly present among the F2 plants used for RNA-Seq, which likely introduced some redundancy among the assembled contigs: more than one contig in ALOMYbase may represent a same unigene and/or several contigs may represent different segments of the same unigene. The precise size of A. myosuroides genome is unknown. Estimations based on DNA cell contents give an expected genome size between 1200 and 4330 Mb [41, 42], i.e., a relatively large genome. Large genomes have a low proportion of transcribed sequences encoding proteins (e.g., 5 to 8 % of the large genomes of fully sequenced legume species genomes [38]). From the range of A. myosuroides genome size estimates, ALOMYbase sequence data would represent 1.6 to 5.7 % coverage of the A. myosuroides genome. Even though redundancy was present in ALOMYbase, this value was deemed acceptable for a transcriptomic resource obtained only from the aerial part of young plants. Five grass species currently have a fully sequenced genome: the weed B. distachyon (270 Mb, 25,532 genes encoding proteins), and the crops O. sativa (rice, 380 Mb, 40,331 genes), S. bicolor (sorghum, 730 Mb, 34,497 genes), H. vulgare (barley, 5100 Mb, c.a. 32,000 genes) and Z. mays (maize, 2100 Mb, 63,540 genes). The genomes of these species contain 95, 106, 47, 6 and 30 genes encoding proteins per Mb genome, respectively. Depending on A. myosuroides genome size estimate, the 32,138 peptides predicted to be encoded by the potentially active contigs in ALOMYbase corresponded to seven to 27 potential genes encoding proteins per Mb genome. Even considering partial transcriptome coverage and occurrence of redundancy, this is in the range of values observed for the grass species with genome sizes similar to that of A. myosuroides (barley and maize).
The N50 size value of the contigs in ALOMYbase was 1240 nucleotides (Table 1), a value higher than those obtained for other recent plant de novo transcriptome assemblies based on Illumina sequence reads [43, 44] or on a combination of Illumina and 454 pyrosequencing reads [45], and similar to that obtained with 454 sequence reads [15]. The average contig size in ALOMYbase (1046 nucleotides) matched the average length of gene coding sequences in grasses (1000 to 1300 nucleotides [46]).
The predicted peptide content of ALOMYbase was compared to those of the five fully sequenced grass genomes. In total, 30.6 % of the protein families identified using OrthoMCL were shared among ALOMYbase and the five grass genomes (Additional file 2: Figure S3). The five grass genomes shared 41.7 % of the protein families identified. These proportions were in agreement with a previous genome-wide study showing that genome peptide contents was largely shared among grass species, including peptide family representation [47]. While ALOMYbase only represents a part of the total A. myosuroides transcriptome, these results suggest a good coverage and a correct representation of the protein families of A. myosuroides genome in ALOMYbase. Considering the number of shared protein families, ALOMYbase content was closer to those of H. vulgare and B. distachyon genomes (Additional file 2: Figure S3). Similarly, rank correlation for the number of predicted peptides assigned to Pfam families was highest between ALOMYbase and B. distachyon genome. The six grass species considered belong to three major subfamilies in the Poaceae: Pooideae (A. myosuroides, H. vulgare and B. distachyon), Panicoideae (Z. mays and S. bicolor) and Ehrhartoideae (O. sativa) [48]. Thus, similarities in protein family representation are consistent with phylogenetic proximity, as already observed [47].
Peptide annotation identified 2017 ALOMYbase contigs potentially encoding peptides assigned to major families involved in NTSR [5, 10]. Peptides annotated as cytochrome P450 or ABC transporters were particularly abundant in ALOMYbase. This could be due to the probable heterozygous status of the plants sequenced. Heterozygosity has been reported to cause redundancy during assembling, especially in fast-evolving gene families like cytochromes P450 [49]. Other possibilities would be gene evolution via duplication and divergence, a process particularly frequent for cytochromes P450 [50], or occurrence of splice variants, which is expected in plants undergoing a stress because alternative splicing is involved in the regulation of stress response [51]. Even with redundancy present in ALOMYbase, our data suggests A. myosuroides genome would be rich in genes encoding proteins potentially involved in NTSR, which could be a reason why NTSR is so widespread and frequent in this species [28].
Particularities in functions such as stress response generally imply particularities in gene expression regulation [47]. A good representation of the transcriptome of the aerial part of A. myosuroides young plants and confirmation of the relevance of RNA-Seq-based expression data using RT-qPCR make ALOMYbase a reliable resource to investigate the transcriptomic response to herbicides inhibiting ALS in A. myosuroides using RNA-Seq data.
ALS inhibitors are among the most broadly and frequently used herbicides. Yet, their effects are still not totally elucidated [30]. Most studies addressing plant response to ALS inhibitors were performed on broadleaved plants, particularly A. thaliana [33, 34] and Pisum sativum [52–54]. Briefly, ALS inhibition is followed by plant growth arrest and the subsequent slow death of treated plants [30]. Application of ALS-inhibiting herbicides causes a rapid induction of specific stress response pathways, including detoxification-related genes [33, 34]. The biosynthesis of branched-chain amino acids (BCAAs: valine, leucine and isoleucine) is interrupted, causing a decrease in the cell contents in free amino-acids, including BCAAs, and a reduction in protein synthesis [30, 55]. This is rapidly followed by an increase in the cellular free amino-acid contents resulting from increased protease-mediated protein degradation and reduced protein synthesis rates [56–58]. Other effects of ALS inhibitors include carbon metabolism impairment leading to an increase in the leaf cell carbohydrate contents and induction of aerobic fermentation [53, 57, 58]. Central energy pathways are also modified in response to oxidative damage [34].
Using a time-course experiment designed as recommended [59] to sample the transcriptome at different times of the day during the early response of A. myosuroides to ALS inhibitors allowed analysis of the part of the transcriptomic response to iodosulfuron + mesosulfuron common to resistant and sensitive plants (Additional file 3: Table S5 and S6). The transcriptional response started at 6HAT with an up-regulation of contigs assigned to thiamin biosynthesis process and oxidative stress response. Thiamin had been proposed to be involved in several abiotic and biotic stress responses, including protection against oxidative stress [60–62]. Oxidative stress associated to ALS inhibitor action had been reported to be transient and moderate, and not a cause for plant death [63], which is consistent with the early and transient response observed here in A. myosuroides plants. From 24HAT on, there was a strong up-regulation of contigs assigned to protein and BCAAs biosynthesis. This is likely a direct consequence of ALS inhibition. Contigs assigned to “Multidrug transport” were up-regulated from 24HAT on, especially multi antimicrobial extrusion proteins (MATE), which is similar to a previous finding of MATE-encoding contigs being up-regulated in A. thaliana in response to three ALS inhibitors [34]. Contigs assigned to “Respiratory gaseous exchange” were also up-regulated. They included seven putative alternative oxidases that are part of the electron transport chain in mitochondria. This may reflect the activation of the alternative respiratory pathway, i.e. aerobic fermentation consecutive to carbohydrate accumulation. This is consistent with previous studies [53, 54, 58, 64, 65]. From 36HAT on, contigs assigned to processes driving gene regulation were up-regulated, while contigs assigned to photosynthesis were down-regulated. This is consistent with plant growth arrest following ALS inhibition. At the latest time-point studied (73HAT), there were few additional changes in the processes enriched in up-regulated contigs. However, there was a drastic increase in the biological processes enriched in down-regulated contigs. In particular, many contigs assigned to various stress response pathways, including oxidative stress response, were down-regulated at 73HAT, as were contigs assigned to oxygen transport and carbon fixation. The effects of ALS inhibitor action on A. myosuroides plants reflected by their transcriptional response were consistent with the literature. They can also be considered to reflect the two first phases of herbicide stress response [10, 66]: the initial shock phase when stress-signalling pathways are triggered (from 6HAT to 24HAT), and the acclimation phase when plant resources are re-oriented towards the establishment of defences (from 24HAT to 73HAT).
Similarities in the response to ALS inhibitors were observed between A. myosuroides and A. thaliana. The expression patterns of the 93 probable ALOMYbase orthologs of A. thaliana genes used as markers for the signature of the response to ALS-inhibiting herbicides were established. ALOMYbase orthologs of A. thaliana marker genes in groups 1a and 1b showed expression patterns remarkably similar to those of their A. thaliana orthologs in both the resistant and the sensitive pools (Additional file 3: Table S3). Similarities in the expression patterns increased with time, with a maximum similarity observed at 73HAT. Our data thus support previous findings that grasses and broadleaved plants share common regulatory mechanisms of gene expression in response to abiotic stresses [67].
From all the foregoing, it is clear that ALOMYbase quantitative and qualitative data are reliable and relevant to study the response to ALS inhibitors in A. myosuroides. Plants resistant to iodosulfuron + mesosulfuron because of NTSR reacted to the stress exerted by the herbicide application, with a part of their transcriptomic response being similar to that of sensitive plants and consistent with the effects of ALS inhibitors described in the literature. This is in contrast to a previous study where, following ALS inhibitor application, no transcriptomic changes were observed in plants resistant because of a mutant herbicide-resistant ALS, while their sensitive counterparts showed extensive transcriptomic changes [33]. Here, A. myosuroides plants resistant because of NTSR were affected by herbicide action in a first step, but ultimately overcame it, as previously proposed [10]. As expected, determinants of the mechanisms allowing NTSR plants to withstand herbicide action are clearly to be sought in the part of herbicide response that is specific to the resistant phenotype.
Overall, there were many more up-regulated contigs than down-regulated contigs in the sensitive plants following herbicide application, while the opposite was observed for the resistant plants (Fig. 4). This could reflect a random response to the herbicide stress in the sensitive plants, while the response in the resistant plants would be more focused and associated with a down-regulation of functions not crucial to overcome ALS inhibitor action. In the sensitive plants, the induction of oxidative stress response seemed more marked, with a specific up-regulation of contigs assigned to this process that could reflect an herbicide stress stronger than in the resistant plants. In the resistant plants, contigs assigned to energy production were up-regulated from 12HAT on, while no significant up-regulation was observed in the sensitive plants. Plant acclimation to stress is an active process that requires extra energy [66]. This observation may thus correspond to an earlier and stronger onset of the phase of acclimation to the herbicide stress in the resistant plants. General analysis of contigs differentially expressed between the resistant and the sensitive plants did not allow identification of biological processes that could directly explain NTSR.
The primary aim of this work was to identify as many as possible of the contigs presents in plants resistant and/or sensitive to ALS inhibitors before herbicide application and at the early stages of herbicide response, in order to establish a transcriptomic resource as comprehensive as possible to allow subsequent studies of A. myosuroides response and NTSR to ALS inhibitors. For this purpose, and because performing RNA-Seq on 42 RNA samples was not affordable, the 42 individual RNA samples studied were sequenced as pooled samples. Transcriptomic data from each pooled sample reflected the average variation in gene expression induced by biological variation among individuals (using three F2 plants per phenotype and per experimental modality) and by the environment (using two clones per F2 plant and per experimental modality). The pooled samples used for RNA-Seq did not contain distinct biological replicates, as is generally recommended to identify contigs differentially expressed among experimental modalities (e.g., [21, 25]). This experimental design was thus not optimal for the purpose of identifying contigs specifically up-regulated in resistant plants that could be candidate NTSR contigs [10], especially because we observed that variation in gene expression could be substantial among the individual samples constituting each pooled sample (Additional file 4: Figure S4; Additional file 7: Figure S7). However, A. myosuroides transcriptomic response to ALS inhibitor action obtained from ALOMYbase expression data was fully consistent with the literature. Despite the limitations inherent to our experimental design, this encouraged us to seek potential NTSR determinants using ALOMYbase expression data. Contigs with a constitutive up-regulation in the resistant plants that was maintained during the early phases of the transcriptional response to iodosulfuron + mesosulfuron were targeted because NTSR is expected to be most effective in avoiding irremediable physiological damage if constitutive, and because contigs that could be used for NTSR detection would be most useful if constitutively differentially expressed in resistant plants (i.e., in the absence of herbicide application).
The five contigs identified herein with a potential link to NTSR were predicted to encode three cytochromes P450 (CYP71A, CYP71B and CYP81D), one peroxidase (Perox2) and one disease resistance protein (DP01). Cytochromes P450 in families CYP71 and CYP81 had been shown to be directly involved in herbicide degradation in crop or model species [68]. In particular, wheat CYP71C6 had been shown to catalyse degradation of a range of ALS inhibitors [69]. CYP71A, CYP71B and CYP81D are thus potentially interesting candidate NTSR contigs in A. myosuroides. The possible roles of Perox2 and DP01 in NTSR are unclear. ALS inhibitors directly cause oxidative stress, but only as a transient side-effect of ALS inhibition [53]. Thus, a direct role of Perox2 in NTSR is dubious. Peroxidases are involved in a variety of plant physiological processes, including abiotic stress response [70]. Perox2 would rather be a NTSR marker, i.e., a contig which expression level is correlated with NTSR but that has no direct role in NTSR [10]. The same applies to DP01, which potentially encodes a peptide with homologies to a protein involved in response to a biotic stress.
The contigs identified herein differed from NTSR genes or candidate transcripts identified in previous studies. Conversely, previously identified NTSR genes or candidate contigs could not be linked to NTSR in the plants studied in our work (Table 3). In particular, the glutathione-S-transferase AmGSTF1 that had been shown to be a major player in A. myosuroides NTSR to herbicides inhibiting ACCase [13, 14] was not identified in our work as a potential candidate contig involved in NTSR to ALS inhibitors. The same applies to the ALOMYbase homologs of other genes or candidate transcripts associated to NTSR to ACCase inhibitors [12, 15] (Table 3). However, these genes or contigs were regulated by herbicide application (Table 3). This suggests that they may be involved in the broad response to iodosulfuron + mesosulfuron, but not directly in NTSR to these herbicides. A. myosuroides response pathways to ACCase inhibitors and to ALS inhibitors, and thus NTSR to each group of herbicides, may be interrelated, as suggested previously [5]. The ALOMYbase homologs of previously identified genes or contigs associated to NTSR to ALS inhibitors [16, 17, 19] were also not involved in NTSR to ALS inhibitors in the A. myosuroides plants studied herein, which is fully consistent with NTSR being underlain by a diversity of mechanisms that vary among species, populations and individuals [5, 6, 9, 10].
In summary, we identified five contigs that are potential NTSR genes or markers. The expression levels of the five contigs varied among the F2 plants analysed: some resistant plants showed a low level of expression for these contigs, and some sensitive plants displayed a high level of expression for some of these contigs. Such variation among individual plants in the expression of genes with a link with NTSR had previously been observed [15, 16, 19, 71]. Principal Component Analysis suggested two groups of NTSR mechanisms could be present in the F2 plants studied, but did not allow separating all resistant plants from all sensitive plants. Our data suggested that the resistant F2 plants studied were not all resistant because of the same NTSR mechanisms, and that not all contigs with a role in NTSR were identified in our study. As all F2 plants studied herein derived from a single parental plant with NTSR to ALS inhibitors, this suggests the occurrence of a set of genes endowing NTSR in the parental plant with NTSR. This is in accordance with previous data showing that NTSR is under polygenic control in A. myosuroides [72]. Further investigation is clearly necessary to confirm the link of the five contigs identified with NTSR, and in particular whether the three putative cytochromes P450 have a degrading activity against iodosulfuron and mesosulfuron.
