Transcriptomic Analysis of Methyl Jasmonate Treatment Reveals Gene Networks Involved in Drought Tolerance at the Early Developmental Stage in Pearl Millet (Pennisetum Glaucum L. Br)


 Water deficit stress at the early stage of development is one of the main factors limiting pearl millet production. A practice to cope with would be to apply hormones to stimulate plant growth and development. Exogenous methyl jasmonate (MeJA) can improve drought tolerance by modulating key pathways, therefore, we assumed that can occur in pearl millet during the early stage of development. To decipher the molecular mechanisms controlling these pathways, RNAseq was conducted in two pearl millet genotypes, drought-sensitive SosatC88 and drought-tolerant Souna3, in response to 200 µM of MeJA. Transcriptomic analysis between the MeJA-treated and non-treated plants revealed 3270 differentially expressed genes (DEGs) in SosatC88 and 127 DEGs in Souna3. Gene ontology (GO) classification assigned most regulated DEGs in SosatC88 to heme binding, oxidation-reduction process, response to oxidative stress and membrane, and in Souna3 to terpene synthase activity, lyase activity, magnesium ion binding, and thylakoid. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis reveals that DEGs in SosatC88 are related to the oxidation-reduction process, the biosynthesis of other secondary metabolites, the signal transduction, and the metabolism of terpenoids, while in Souna3, DEGs are related to the metabolism of terpenoids and the energy metabolism. Two genes encoding a diterpenoid biosynthesis-related and a Glutathione S transferase T3 were contra-regulated between SosatC88 and Souna3. Additionally, five random genes differentially expressed by RNAseq were validated using qPCR. These insights into the molecular mechanisms of pearl millet genotype tolerance at the early stage of development contribute to the understanding of the role of hormones in adaptation to drought-prone environments.


Introduction
Worldwide, water de cit contributes to crop production losses, particularly in the Sub-Saharan Africa countries, where agriculture depends on rainfall. In addition, population growth increases food demand while arable land is continuously shrinking. Meanwhile, climate forecasts are predicted to be extremes by 2050 1 in Africa, therefore, there is a need to breed high-yielding drought-tolerant genotypes. One of the major staple food and drought-tolerant plant is pearl millet (Pennisetum glaucum L. Br), ranked fourth among the cultivated tropical cereals. Since its domestication from West Africa 4,900 years ago 2 , pearl millet constitutes a staple food for more than 100 million people, with 90% grown by smallholder farmers 3 , hence contributing to nutrition and food security. Indeed, its grain is gluten-free, rich in vitamins, micronutrients (iron and zinc), and protein (8-19%) 4 . With a genome size estimated to 1.79 gigabytes, recent studies led to the annotation of 38,579 genes, identi cation of speci c marker-traits linked to key agronomic traits, such as grain yield, stem and leaf biomass as well as tremendous adaptive traits for heat and drought 5 . As a cross-pollinated C4 plant (2n = 2x = 14), pearl millet endows with a large genetic diversity and grows in some of the hottest and driest areas across the arid and semi-arid regions. Its vulnerability to drought is genotype dependent 6 . Drought-related stresses at early stage of pearl millet plant development are limiting factors to its productivity. Water de cit for example induces a limitation of photosynthesis activities due to reduced leaf expansion, impaired photosynthetic machinery, premature leaf senescence 7,8 , stomatal closure 9 , and the activity of Calvin's cycle enzymes such as Rubisco and phosphoenolpyruvate carboxylase 10 .
Reactive oxygen species (ROS), produced under water de cit, target various organelles, including chloroplasts, mitochondria, and peroxisomes, resulting in cell membrane instability, senescence, or plant death 11 . Other physiological responses to water de cit include decreasing stem height and diameter, leaf number, leaf area ratio, dry matter, shoot/root weight ratio, net CO 2 assimilation and chlorophyll uorescence, photosynthetic activity, altered cell wall elasticity and generation of toxic metabolites 12,13 .
Another important factor leading to plant tolerance to water de cit is hormonal regulation. The expression of several transcription factors and their target genes involved in mediating phytohormones abscisic acid (ABA) under water de cit conditions, key components of perception and signaling, are induced to attenuate the negative impacts due to water de cit 12,14 . Exogenous treatment of abscisic acid 15 and salicylic acid 16 improves tolerance under water de cit in wheat (Triticum aestivum). The application of jasmonates induced improved drought tolerance in barley (Hordeum vulgare), 17 , soybean (Glycine max) 18 , and sugar beet (Beta vulgaris) 19 . A treatment of Methyl jasmonate (MeJA) can enhance tolerance to dehydration of plants subjected to osmotic stress caused by polyethylene glycol 20 and water retention 21 . MeJA governs many aspects of plant development, including seed germination, root growth, owering, fruit ripening, senescence 22,23 and abiotic stresses tolerance 24 . This hormone mitigates the in uence of water de cit by promoting stomatal closure induced by H 2 O 2 generation 25 , acting principally on the expression of genes involved in signaling processes.
Transcriptomic studies in pearl millet have been used to reveal responsive genes under drought, heat, or salinity conditions 26,27 . However, the mechanism of pearl millet in response to MeJA treatment has a negligible focus. The present study investigates genes that are differentially expressed after an exogenous MeJA treatment and with a mimic of water de cit by withholding water supply to decipher the molecular basis modulating the response of pearl millet under water scarcity circumstances. By spraying pearl millet plants at the early stage of development with MeJA, the assumption that a differential gene expression could activate or repress genes underlying key metabolic and biosynthesis pathways that lead to coping with water de cit-related stress was tested. Comparative transcriptome analysis of two pearl millet genotypes contrasting by their level of tolerance to water de cit at the early stage identi ed sets of differentially expressed genes (DEGs) between non-treated and MeJA-treated plants. The regulation of DEGs modulated by MeJA and new genetic insights of adaptive mechanisms of pearl millet in response to water de cit stress are discussed.

Results
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 signi cant 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 signi cantly change leaves number in SosatC88 at T1, i.e ve 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 signi cantly 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 signi cant 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 de cit and remained lower in MeJA-treated plants, however, this change was only signi cant 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 gures, Fig. S1), indicating that MeJA penetrated pearl millet plants and did change gene expression. A signi cant positive fold change PgNCED expression pro le was observed in both droughtsensitive 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 gures, 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 ltration 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 nal number of 367,493.828 high-quality reads (passing lter 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 le, Table S1).
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 identi ed above. GO classi cation 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 classi cation 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 signi cantly enriched (Fig. 5a). Between MeJA-treated SosatC88 and Souna3, none of the categories are signi cantly 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 signi cantly 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_MeJAtreated 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 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, stressresponsive 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 signi cant correlation (r 2 = 0.777, p = 0.02) was observed between qPCR and RNA-Seq data that validate the differential expression of genes (Fig. 8c).

Discussion
Water de cit is considered one of the most threatening abiotic stresses in agriculture in the era of climate change 28 . Water de cit resulting from these drought periods affects different aspects of plant growth and development, leading in some cases to a considerable crop production loss. These negative impacts can be alleviated by a wide range of plant hormonal treatments. MeJA, which is known to be involved in stress resistance mechanisms in plants, may have a key role to play in providing useful solutions. MeJA signi cantly reduced growth and leaf expansion in soybeans (Glycine max), hampered seedling diameter in scots pine (Pinus sylvestris) [29][30][31] , and inhibited both root and shoot growth in rice (Oryva sativa) 32 . The treatment of MeJA 200µM at the early stage of pearl millet development did not induce signi cant visible morphological change between genotypes. A reduction in chlorophyll content suggests a loss of photosynthetic pigments or activity. These are in accordance with previous studies 30,33,34 , that observed a decline in photosynthesis of soybean (Glycine max), barley (Hordeum vulgare) leaves and Arabidopsis thaliana. The loss of photosynthetic pigments would decrease the amount of energy absorbed by the photosynthetic system, thus reducing the energy required for metabolism, which consequently impacts plant growth and development. The loss of chlorophyll content could also be attributed to leaf senescence 35 . Furthermore, it has been shown that jasmonates attenuate the in uence of drought by favoring the closure of stomata 25,36 . Our ndings suggest that treatment of MeJA could cause the stomatal closure, leading to the inhibition of CO 2 absorption and the reduction of photosynthetic activity 10,37,38 .
Despite non-signi cant morphological changes between genotypes, MeJA treatment was effective in modulating endogenous stress-related genes expression. It is known that gene encoding Lox2 is highly induced by jasmonates (JA and MeJA) or in response to biotic stress 39 and NCED is the key enzyme in the ABA biosynthetic pathway in plants 40,41 . The overexpression of pearl millet Lox2 (PgLox2) is induced by the jasmonates produced by the plants in response and in response to MeJA treatment. The expression pro le of pearl millet NCED (PgNCED) in MeJA-treated and non-treated plants suggests that water de ciency and MeJA treatment induced an increase in ABA biosynthesis, which is most important with external MeJA input. These results suggest that the decrease in plant growth is caused by increased ABA and MeJA accumulation.
To gain insights in the molecular basis of pearl millet responses to MeJA treatment, a transcriptome analysis was carried out at the early stage of plant development. Transcriptomic data were generated to identify the molecular basis at the genome-wide level underlying the response of pearl millet to MeJA treatment at the early developmental stage. Under non-treated conditions, distribution of DEGs in droughtsensitive SosatC88 and drought-tolerant Souna3 in extracellular region term, pathogenesis term and transferase activity, transferring acyl group term according to GO analysis suggested that waterde ciency involves genes related to these terms. In contrast, when MeJA was applied there was no signi cant difference in the transcriptomic response between SosatC88 and Souna3. This suggests that MeJA stimulates in SosatC88 pathways and biological processes to better adapt to water de ciency as in Souna3. In addition, 12 genes were found to be differentially expressed in the opposite manner depending on whether we are in MeJA-treated or non-treated conditions. Among these genes, two are upregulated by the MeJA treatment. These two genes are Pgl_GLEAN_10004431 encoding for an unknown protein and Pgl_GLEAN_10018402 encoding for an uncharacterized protein family UPF0136 transmembrane. The additional ten genes were downregulated by the MeJA treatment. They are mostly Tyrosine-protein kinase (Pgl_GLEAN_10005256), Multicopper oxidase, type 2 (Pgl_GLEAN_10001778), Peptidase (Pgl_GLEAN_10016112), Phospholipid/glycerol acyltransferase (Pgl_GLEAN_10012582), Glycoside hydrolase (Pgl_GLEAN_10004142) and uncharacterized genes. About the uncharacterized genes, they could likely represent key regulators of the variation in owering time or drought tolerance among different genotypes. Further investigation is required for understanding the role in pearl millet adaptation to stress.
Data reported here indicate that the MeJA triggers higher transcriptional responses in SosatC88 to better adapt and cope with water de ciency in opposite to Souna3. The striking difference in DEGs between SosatC88 and Souna3 provides information on the amplitude of the effect of the MeJA treatment on the transcriptome of SosatC88. Consistent with data reported here, other transcriptome studies report a higher number of DEGs under water stress-sensitive variety in wheat (Triticum aestivum) 42 and sorghum (Sorghum bicolor) 43 . MeJA induced in both genotype enrichment of GO terms related to terpene synthetase activity and lyase activity. One way to protect against oxidative stress is the biosynthesis of volatile terpenoids, which is thought to quench ROS 44 . Lyases are known to play a role in antioxidant activity in cases of biotic or abiotic stress 45 . We theorize that these groups of genes probably constitute the characteristic basis for the adaptation to water de cit conferred by the treatment of MeJA in both genotypes SosatC88 and Souna3. On the other hand, the results of the GO enrichment analysis of upregulated gene groups provide evidence of striking antioxidant enzyme activities. Several important functional gene groups related to water de ciency stress including oxidation-reduction process, heme binding, electron transfer activity, response to oxidative stress, membrane, peroxidase activity, among others, are signi cantly enriched in SosatC88 treated with MeJA. Peroxidases are involved in the defense against abiotic stress through their role in the scavenging of ROS 46 . They are also involved in cell wall formation, confer cell wall rigidity and reducing leaf expansion 29 . It has also been shown in several studies that MeJA treatment increases peroxidase activity 29,34,47 . A cross-species meta-analysis of progressive drought stress identi es the heme-binding family protein as key drought adaptive genes conserved 48 . These results show that MeJA improves SosatC88 adaptation, drought-sensitive genotype, to water de ciency by increasing and diversifying its antioxidant response mechanisms to cope under these unfavorable conditions.

Conclusion
A successful treatment of exogenous MeJA in two genotypes with different response to drought did not alter plant morphology, but modulated gene expression. RNAseq analysis revealed a transcriptional process that mobilizes a large number of DEGs related to cell wall protection and detoxi cation, implying different mechanisms involved in response to drought at the early stage of the pearl millet developmental cycle. Five genes differentially expressed are genuine candidates as molecular markers in functional genomics or in the breeding of pearl millet. Overall, data give an indication that, in pearl millet, exogenous MeJA treatment may play a dual role by simulating abiotic stress defenses to protect plants under water de cit and contribute to crop production in drought-prone environments.

Plant material and growth Conditions
Seeds of two pearl millet open-pollinated varieties, SosatC88 and Souna3, were used. SosatC88 was developed by recombining 19 S1 progeny selected at Cinzana (Mali) in 1988 from a cross between the local Souna and Sanio genotypes. It was considered a water de cit-sensitive check. Souna3 is a synthetic variety resulting from the recombination of eight lines from the PC 28 and PC 32 populations (106-7, 108-4, 113-3, 115-4, 134-5, 142-4, 143-4, 148-3). Considered as a water de cit tolerant genotype, Souna3 reaches maturity between 85-95 days and can be grown in areas with moderate rainfall ranging between 400 and 750 mm was. All plant material used in this study came from ISRA genebank and their uses comply with relevant institutional, national, and international guidelines and legislation. Both genotypes were sown in 10 L buckets containing a well-characterized soil from Bambey (Senegal) 49 . In each bucket, 4 plants of the same variety were grown. These buckets were completely randomized with 6 replicates per treatment.

Page 10/19
The plants were watered three times a week by pouring an equivalent volume of water on each soil bucket until 17 days after sowing (DAS). Then, water was withdrawn and plant leaves were sprayed with a solution containing 200 µM MeJA (Methyl jasmonate 95%, CAS Number: 39924-52-2 (Methyl 3-oxo-2-(2pentenyl) cyclopentaneacetate from Sigma-Aldrich) dissolved in 0.1% tween 20 and with 0.1% tween 20 for the MeJA-treated and the non-treated, respectively, every day for 10 days. Three days after the last treatment, the aerial parts of the MeJA-treated and non-treated plants were collected and immediately frozen in liquid nitrogen and stored at -80° C for further RNA extraction.

Development and Chlorophyll density parameters measurement
Leaf number and plant height were assessed at T0, T1 and T2. T0 corresponds to the rst day just before treatment (both non-treated and MeJA-treated plants), T1 and T2 correspond to 5 and 10 days after treatment, respectively. The chlorophyll density level was measured at T1 and T2 using a SPAD 502 PLUS CHLOROPHYLL-METER. Four plants from the same pot were used for measurement of leaves number, plant height and chlorophyll density level, which resulted in six biological replicates for each genotype in each sampling date. A signi cant difference in mean between the samples as obtained by Student's t-test at P < 0.05.

RNA extraction
Leaves from four plants of each genotype were pooled in one sample. Sampling was performed on three biological replicates (MeJA-treated and non-treated plants). The samples were ground in liquid nitrogen and then, 200 mg of powder was used for total RNA isolation with TRI Reagent (Sigma-Aldrich) according to the instructions of the manufacturer. RNA was then DNase treated (Qiagen, Germany) and puri ed using Qiagen RNeasy Minikit. The quality of the RNA was evaluated on 1.2% agarose gel (Supplemental gures, Fig. S2) and electrophoresis was carried out at 100 Volts for 30 min and the quanti cation was performed using a Nano Drop Lite and the Agilent 2100 bioanalyzer.

Gene markers expression level assay for MeJA treatment validation
To validate the results observed at the end of the MeJA treatment, we identi ed genes that are involved in response to plant biotic and abiotic stresses 23 . Among these, we used Lipoxygenase 2 (Lox2) involved in the biosynthesis of MeJA and 9-cis-epoxycarotenoid (NCED) involved in abscisic acid (ABA) biosynthesis. Both are associated in plant response to abiotic/biotic stresses. Lox2 encodes the enzyme that catalyzes an important step in jasmonic acid (JA) biosynthesis from linolenic acid derived from the membrane 50,51 while NCED is the key enzyme in ABA biosynthesis involved in plant response to water de cit 41 .
Primer pairs for PgLox2 and PgNCED were designed using Primer3Plus online tool 52 according to the target gene-speci c sequence design, the option of qPCR setting Primer3Plus online tool platform was used. The primer sequences of transcripts are listed in Supplemental le, Table S2. qPCR assays Five hundred (500) ng of the total RNA was reverse transcribed using the Goscript Reverse Transcriptase system (Promega) according to the manufacturer's instructions. The qPCR reaction was performed in a total volume of 10 µl containing 5 µl of GoTaq® qPCR Master Mix (2X), 0.25 µM of each primer, 1 µl (100 ng) of diluted cDNA and free nuclease water and then, loaded in a StepOnePlus™ Real-Time PCR Systems. The program used for qPCR contains three steps: the rst step is the initial activation at 95°C for 2 min ; the next step is the PCR ampli cation performed up to 40 cycles with 95°C for 15 s, 58°C for 60 s and the third step is a melt curve analysis ramping from 58°C to 95°C.
The relative expression levels of transcripts under the different conditions were normalized based on the expression of an endogenous control gene, ubiquitin 5 (UBQ5) 53 using the 2 −∆∆Ct method 54 . The qbase+ software Version: 3.2 55 was used for analysis.

RNAseq
Library construction and sequencing For library construction, RNAs were puri ed on RNAClean XP beads and then analyzed by capillary electrophoresis (Fragment Analyser). Highly pure mRNA was isolated from the total RNA using oligo (dT) beads. The Illumina TruSeq RNA Library Prep Kit v2 was used to synthesize the second strand cDNAs library. Sequencing was performed on an Illumina HiSeq 2500 using the Sequence By Synthesis (SBS) technique using the TruSeq Rapid SBS kit (Illumina). Sequencing results were obtained as single-end reads (50 bp each) in the FASTQ format.

Differential gene expression analysis
Clean reads were aligned to the reference genome 5 and counted with Bioconductor R package Rsubread 56 . The read counts were used to performed the Differential Expression Genes analysis with edgeR package of Bioconductor 57 . Pairwise comparisons were performed between four (4) datasets: SosatC88 vs Souna3(non-treated)", SosatC88 vs Souna3 (MeJA-treated)", SosatC88 MeJA-treated vs SosatC88 non-treated" and Souna3 MeJA-treated vs Souna3 non-treated". False Discovery Rate (FDR) < 0.05 and log|Fold change| > 2 were considered as the conditions to state the genes as differentially expressed.

Functional annotation of transcripts and genes enrichment analysis
The sequence similarity search was carried out against the National Center for Biotechnology Information (NCBI) non-redundant (nr) protein database and the Swiss-Prot protein database using the BLASTx algorithm from Blast2GO 58 specifying E-values < 10 −5 . Gene Ontology (GO) categorization was done with Blast2GO. From the annotated GO le, we classi ed and plotted the DEGs according to the o cial classi cation (Molecular Function, Biological Process and Cellular Component), using WEGO 2.0 59,60 and functional enrichment analysis was also performed using Goseq, a bioconductor R package implemented in galaxy 61 using Wallenius method.
DEGs were then mapped on the biological pathways of the web-based Kyoto Encyclopedia of Genes and Genomes (KEGG) web-based annotation server (KASS) by running Blastx against KEGG GENES (Kyoto Encyclopedia) https://www.genome.jp/kegg/kaas/. We completed the KEGG annotation using the KEGG result les of millet genome information uploaded to GIGAdb 62 . The result contains KO assignments (KEGG Orthology) and KEGG pathways were automatically generated. KEGG enrichment analysis was also performed based on the same criteria described above for GO.

Declarations
The authors declare no competing interests.

Data availability
All data provided in the article and Supplementary les are available from the corresponding author on reasonable request.     The top 10 over-represented GO terms under non-treated condition between SosatC88 and Souna3. b) The top 10 over-represented GO terms under MeJA-treated condition between SosatC88 and Souna3. The Adj p-value is the corrected p-value ranging from 0 to 1.   Validation of RNASeq data with qPCR. a and b) Expression of randomly selected genes was examined by qPCR analysis. For each gene, fold changes were calculated by ΔΔCt method and log2Fold change were compared between qPCR and RNAseq. c) Correlation between RNAseq and qPCR data based on log2fold change of the ve selected genes: y=0.31x -0.8.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.