MeJA effects at early-stage development of pearl millet
Plants of two genotypes contrasting for the response to drought (sensitive SosatC88 and tolerant Souna3) showed no significant developmental changes for leaves number and plant height, p > 0.05 at 17 days after sowing, i.e T0. The aerial part of these plants was then daily sprayed with MeJA (200 µM) for 10 days over T0 and we observed that the non-treated plants of both SosatC88 and Souna3 appeared withered compared to those treated (Fig. 1a,b). Plant height, leaves number and chlorophyll content were measured in both genotypes. MeJA did not significantly change leaves number in SosatC88 at T1, i.e five days after treatment (p=0.72), and T2, i.e 10 days after treatment (p=0.052). However, in Souna3 at T1 MeJA-treated plants have significantly more leaves (p=0.049) than non-treated plants (Fig. 1c,d). In both genotypes, the height of MeJA-treated plants stayed lower than that of non-treated plants until the end of the treatment with a significant difference observed (p=0.03) in SosatC88 at T2 (Fig. 1e,f) in treated and non-treated plants. In both genotypes, the chlorophyll content decreased during the 10 days of water deficit and remained lower in MeJA-treated plants, however, this change was only significant in Souna3 at T1 (p=0.03) (Fig. 1g,h).
MeJA treatment in modulating gene expression
To evaluate whether the MeJA treatment in both SosatC88 and Souna3 mimicked drought tolerance response at the gene expression level, we conducted qPCR using RNA extracted from the aerial part of plants sprayed or not for 10 days with MeJA. The relative expression levels of pearl millet lipoxygenase 2 (PgLox2) and of 9-cis-epoxycarotenoid (PgNCED) were determined. The two genes were previously used as endogenously regulated after a MeJA treatment or in response to drought-induced stress, respectively. The results showed that PgLox2 was over-expressed in MeJA-treated plants compared to the non-treated (Supplemental figures, Fig. S1), indicating that MeJA penetrated pearl millet plants and did change gene expression. A significant positive fold change PgNCED expression profile was observed in both drought-sensitive SosatC88 and drought-tolerant Souna3, indicating that the MeJA treatment modules the response to drought-induced stress in pearl millet by modulating gene expression (Supplemental figures, Fig. S1).
RNA-seq assembly and analysis after MeJA treatment
To further explore this contrasting behavior between genotypes, we analyzed the transcriptomic responses to MeJA treatment on drought-sensitive SosatC88 and drought-tolerant Souna3. A total of 408,935.888 raw reads were generated, with an average of 34,077.991 reads per sample. After filtration and quality check, 379,364.208 high-quality reads were obtained with an average of 30,624.486 reads per sample. On average, 92.77% of total data passed a Phred score of ≥ Q30 bases for each sample. A final number of 367,493.828 high-quality reads (passing filter and correctly assigned to a sample after demultiplexing) were then aligned onto the pearl millet reference genome and 89.29 – 90.81% reads from twelve samples were uniquely mapped (Supplemental file, Table S1).
Pairwise comparison (Fig. 2a-c) between non-treated SosatC88 and Souna3 plants revealed a total of 517 genes upregulated and 149 genes downregulated (dataset 1). Between MeJA-treated SosatC88 and Souna3, 265 genes are upregulated and 473 genes are downregulated (dataset 2). When comparing within the genotype, pairwise analysis shows 3270 DEGs (2553 upregulated and 717 downregulated) for SosatC88 (dataset 3) and 127 DEGs (71 upregulated and 56 downregulated) for Souna3 (dataset 4). Comparison from dataset 1 and dataset 2 revealed that 70 and 40 genes are upregulated and downregulated, respectively, however, 12 genes were contra-regulated (Fig. 2d). Comparison from dataset 3 and dataset 4 showed 49 and 25 genes upregulated and downregulated, respectively. However, two genes encoding a diterpenoid biosynthesis-related and a Glutathione S transferase T3, respectively, were contra-regulated between SosatC88 and Souna3 (Fig. 2e).
DEGs functional analysis after MeJA treatment
Gene ontology (GO) assignments of biological process (BP), cellular component (CC) and molecular function (MF) were used to classify the differentially expressed genes and their expected functions from the datasets identified above. GO classification of pearl millet gene lists from dataset 1, dataset 2, dataset 3 and dataset 4 were also carried out. From dataset 1, a total of 45 GO terms were assigned (13 CC, 10 MF and 22 BP) (Fig. 3a), 40 GO terms for dataset 2 (13 CC, 9 MF and 18 BP) (Fig. 3b), 48 GO terms for dataset 3 (14 CC, 10 MF and 24 BP) (Fig. 4a), and 30 GO terms for dataset 4 (10 CC, 4 MF and 16 BP) (Fig. 4b). This classification reveals that in non-treated conditions, differences between the two genotypes are particularly in cell (141 genes) cell part (140 genes), membrane (122 genes), catalytic activity (273 genes), binding (235 genes), metabolic process (252 genes) and cellular process (176 genes). In MeJA-treated conditions, the most differences appeared in cell (179 genes), cell part (177 genes), membrane (191 genes), membrane part (149 genes), catalytic activity (305 genes), binding (294 genes), cellular process (223 genes), and metabolic process (273 genes) indicating that some important cellular processes and metabolic activities occurred in the leaves of water-stressed millet in response to MeJA treatment.
According to the GO ontology, MeJA treatment modulated several genes involved in cell part, cell, organelle, protein-containing complex, organelle part, membrane, membrane part, extracellular region, catalytic activity, binding, transporter activity, localization, biological regulation, response to stimulus, metabolic process, cellular process, regulation of the biological process, signaling, cellular component organization or biogenesis. In drought-sensitive SosatC88, the DEGs are greater than in drought-tolerant Souna3 (Fig. 4a,b). However, MeJA treatment during non-watered conditions has a particular impact on the genes involved in the structure of the membrane and the binding genes.
GO terms enrichment analysis
Goseq analysis was used to classify the functions of DEGs. Between both non-treated SosatC88 and Souna3, with related terms ‛extracellular region’ (cellular component category), ‛pathogenesis’ (biological process category) and ‛transferase activity, transferring acyl group’ (molecular function category) were significantly enriched (Fig. 5a). Between MeJA-treated SosatC88 and Souna3, none of the categories are significantly enriched; but molecular function category including ‛magnesium ion binding’, ‛lysase activity’, ‛terpene synthase activity’, ‛DNA binding’, ‛nucleotide binding’, ‛polysaccharide binding’, ‛calcium binding’, ‛oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor’ and cellular component including ‛nucleus’ and ‛nucleosome’ were the top over-represented categories (Fig. 5b).
Between SosatC88 MeJA-treated and non-treated, molecular function category including ‛heme binding’, ‛iron ion binding’, ‛terpene synthase activity’, ‛lyase activity’ ‛oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen’, ‛electron transfer activity’, ‛peroxidase activity’ and biological process category including ‛oxidation-reduction process’, ‛photosynthesis, light harvesting’, ‛response to oxidative stress’ were the most highly enriched (Fig. 6a). On the other hand, molecular function category, ‛terpene synthase activity’, ‛lyase activity’, ‛magnesium ion binding’; cellular component category, ‛thylakoid’ and biological process ‛photosynthetic electron transport chain’ were significantly enriched in Souna3_MeJA-treated vs non-treated. The GO molecular function category, ‛terpene synthetase activity’ and ‛lysase activity’ were also highly enriched in both SosatC88_MeJA-treated vs non-treated and Souna3_MeJA-treated vs non-treated (Fig. 6b).
Pathways analysis
Based on the KEGG database, several signaling pathways and genes involved in the response of both millet varieties to MeJA treatment in water-stressed conditions were discovered. Comparing dataset 1 and dataset 2, an increase in carbon metabolism, Pyruvate metabolism, Oxidative phosphorylation was observed (Fig. 7a,b and Supplementary file). Propanoate metabolism, Inositol phosphate metabolism, Carbon fixation in photosynthetic organisms, Biosynthesis of unsaturated fatty acids are induced by MeJA treatment. The most representative pathways of the overexpressed genes in dataset 3 and dataset 4 were those associated with metabolism (KO 01100), biosynthesis of secondary metabolites (KO 01110), antibiotic biosynthesis (KO 01130), carbon metabolism (KO 01200), photosynthesis (KO 00195) and plant hormone signal transduction (KO 04075) (Fig. 7c,d and Supplementary file). Some enriched pathways such as oxidative phosphorylation (KO 00190), carbon fixation in photosynthetic organisms (KO 00710), phenylalanine metabolism (KO 00360), phenylalanine tyrosine tryptophan metabolism (KO 00400), phenylpropanoid (KO 00940), ubiquitin-mediated proteolysis (KO 04120), MAPK signaling (KO 04016) and circadian rhythm (KO 04712) are only represented in SosatC88, indicating that these pathways were recruited by the MeJA hormone to enhance SosatC88 adaptively under water deficit.
Validation of gene expression using qPCR
Five genes differentially expressed under or not treatment of MeJA of SosatC88 and Souna3 were selected based on log2 fold change and adjusted p-value parameters (FDR < 0.05) and their expression was validated by qPCR. These DEGs encode for: a small Auxin Up-Regulated RNA 50 (SAUR50), a peroxidase A2 (PA2), a myeloblastosis 30 (MYB30), a receptor kinase-like protein Xa21 and a glutathione S transferase T3 (GstT3). These DEGs included plant hormone signal transduction related genes, stress-responsive genes, peroxidase activity-related genes. Differential expression of transcripts was observed in SosatC88 and Souna3 between non-treated and MeJA-treated plants in concordance with the RNAseq data (Fig. 8a,b). The results also validated the contra-regulation revealed in RNAseq of PgGstT3 by MeJA treatment in SosatC88 and Souna3. A Pearson correlation was performed to test whether the qPCR results correlated with the RNAseq results. The ratio of expression levels obtained between non-treated and treated plants by qPCR was compared to the ratio of expression measured by RNA-Seq. A significant correlation (r2 = 0.777, p = 0.02) was observed between qPCR and RNA-Seq data that validate the differential expression of genes (Fig. 8c).