Transcriptomic studies reveal a key metabolic pathway contributing to a well-maintained photosynthetic system under drought stress in foxtail millet (Setaria italica L.)

Materials and experimental design

The materials used in this study were E1 (maternal line), H1 (paternal line) and their F₁ hybrid M79, a drought resistant variety. H1 is a drought-tolerant cultivar released from the Shanxi Academy of Agricultural Sciences, and has been sporadically grown in arid and barren areas during the past decade. E1 developed from the same institute has been an important parental line for many varieties currently grown in China.

Millet seeds were surface sterilized with 0.5% NaClO, washed three times with ddH₂O and sown in pots (diameter: 8 cm, height: 10 cm, 15 seeds per pot) filled with peat and nutrient soil (1:1). After the pots were kept in a growth chamber (light/darkness: 16h/8 h, temperature: 30 °C/22 °C) for seven days, five healthy plants were maintained in each pot by removing extra plants. After three weeks, the plants at six-leaf stage were treated with drought stress. Each genotype with 100 plants were divided into the control group (M79_CK, E1_CK, and H1_CK) with regular watering, and the drought group (M79_DR, E1_DR, and H1_DR) without watering. For the latter group, the soil gravimetric water content was monitored using a soil moisture analyzer TDR300 (Spectrum, Aurora, IL, USA). At the 5th, 6th, and 7th days after stopping watering, the soil gravimetric water content dropped to 26.1%, 20.3%, and 14.6% of field capacity, respectively.

Samples were collected on the 6th day when the soil water content was about 20% (Tang et al., 2017). At this stage, the water potential in M79 leaves was significantly higher than those of the parental lines (Fig. 1B). Plants from the control group (soil water content about 55%) were also sampled at the same stage. Samples were immediately wrapped in foil and frozen in liquid nitrogen, and then kept in a −80 °C freezer. The top second leaf was used as the test materials for all treatments with the upper half for transcriptome sequencing and the lower half for measuring physiological indicators including catalase (CAT) and relative electrolyte leakage (REL). All samples were biologically duplicated three times. The top second leaves from different plants of the same treatments were used for measuring leaf water potential (LWP), photosynthetic rate (A), transpiration rate (E), maximum energy conversion efficiency in PSII centers (Fv/Fm), and quantum yield of PSII electron transport (ΦPSII). The back of each tested leaf was labeled before measuring LWP, A, E, Fv/Fm and ΦPSII to facilitate the next measurement. Each measurement was repeated five times using fully expanded, uninjured leaves.

Measurement of drought-related physiological changes

Plasma membrane damage of millet leaves was determined as previously reported (Cao et al., 2007), and REL was used to measure the extent of damage to the plasma membrane. CAT quantification was performed following Bonnecarrère et al. (2011). LWP was measured using a Psypro plant water potential meter (WESCOR, Logan, UT, USA). A, E, Fv/Fm, ΦPSII were measured using a LI-6800 portable photosynthesis system (LI-COR, Lincoln, NE, USA) as described by Lowry et al. (2015). Light intensity, CO₂ concentration, and air flow rate were set to 800 µmol m⁻² s⁻¹, 400 µmol mol⁻¹, and 500 µmol s⁻¹, respectively. Measurements were carried out from 8:30 - 11:30 am on each day. WUE (µmol mmol⁻¹) was calculated by: WUE = A (µmol m⁻² s⁻¹)/E (mol m⁻² s⁻¹)/1000.

RNA extraction, cDNA library construction, and transcriptome sequencing

RNA samples were prepared from 18 harvests (2 treatments × 3 genotypes × 3 biological replicates) using Trizol^® reagent (Invitrogen, Waltham, MA, USA) for subsequent RNA-seq analysis. Quality and concentration of RNA were determined by agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). A more accurate RNA quantification was performed by using an Agilent 2100 Bioanalyzer RNA Nano chip (Agilent Technologies, Waldbronn, Germany).

Messenger RNA from each sample was enriched using mRNA Capture Beads and then fragmented at high temperature in the presence of metal ions. Using mRNA as a template, the first strand cDNA was synthesized using random hexamers, followed by synthesis of the second strand cDNA, and then purification of the double-stranded cDNA using VAHTS™ DNA Clean Beads. After end repair and A-tailing, purified double-strand cDNA was ligated to the sequencing adapter and sorted using VAHTS™ DNA Clean Beads to get 300–400 bp fragments. Finally, PCR amplification was performed, and PCR products were purified with VAHTS™ DNA Clean Beads to generate the final libraries. Library concentrations were assayed using Qubit 3.0, and the library inserts were subsequently examined using an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and further quantified using the ABI Step One Plus Real-Time PCR system. Finally, the libraries were pooled and sequenced on a HiSeq X Ten (Illumina, San Diego, CA, USA) platform using PE150 mode by Nanjing Vazyme Biotech Company, Ltd.

RNA-seq data analysis

Clean reads extracted from raw reads using Tophat2 (v2.0.13) were compared with the reference genome Setaria italica V2.2 (phytozome.jgi.doe.gov) to get mapped reads (Kim et al., 2013). Based on the available data, we also performed analyses on gene saturation, homogeneity of sequencing, the proportions of mapped reads in genomic exons, introns, and intergenic regions, and correlation analysis between replicates.

Gene expression analysis was performed using Cufflinks v2.2.1. Transcriptome reads aligned to the reference genome were quantified and normalized to fragments per kilobase of transcript per million fragments mapped (FPKM), differences between drought-treated, and the control FPKM values were compared using the software Cuffdiff v2.2.1 (Trapnell et al., 2010). The thresholds of DEGs were set as FDR ≤ 0.05 and |log₂ FoldChange|≥ 1.

Functional annotation, pathway analysis, clustered heat map, and co-regulation network analysis

All DEGs were mapped to terms in the GO database (http://www.geneontology.org/), and the number of genes per term was calculated. Based on the GO:: TermFinder, the GO enrichment of the DEGs was performed using a hypergeometric test with a corrected FDR < 0.05 as a threshold (Boyle et al., 2004). A biological pathway analysis of DEGs was performed using KEGG (http://www.genome.jp/kegg/), and significance was calculated by hypergeometric distribution with a corrected FDR < 0.05 as the threshold (Kanehisa et al., 2008).

To generate a clustered heat map, expression data were converted using the formula log₂ (FPKM + 1), and the map was drawn using the heatmap 2 function in the R/Bioconductor package gplots (Warnes, 2016).

Co-regulation network analysis was conducted by using Cytoscape (v3.4.0) to plot the co-regulation network with the Pearson correlation coefficient setting |PCC| ≥ 0.93 (Saito et al., 2012).

qRT-PCR

To verify RNA-seq data, DEGs were confirmed by qRT-PCR following Livak & Schmittgen (2001). Primers were designed (Table S2) based on gene sequences from Setaria italica V2.2 (phytozome.jgi.doe.gov) using Primer 3 (http://frodo.wi.mit.edu/). Quantitative PCR was performed using a SYBR^® Green PCR Master Mix Kit (Applied Biosystems, Foster City, CA, USA) and an ABI7900 system. The 2^−ΔΔCT method was used to calculate relative gene expression. Correlation between RNA-seq and qRT-PCR was analyzed using SPSS 22.0 software (IBM, USA).

Phenotypes of M79 and its parental lines under drought conditions

After drought treatment for five days, seedlings of M79 and its parental lines showed no obvious differences. However, after 6 or 7 days of drought stress, leaves in most E1 and H1 plants exhibited curling and withering, whereas M79 appeared normal (Fig. 1A). Physiological analysis of those plants before and after drought treatment showed that LWP decreased sharply after drought stress. On the 5th day of drought treatment, the water potentials of E1 and H1 significantly declined (P < 0.01 and 0.05, respectively) in comparison to M79 (Fig. 1B). We used REL to analyze the amounts of damage after drought. From the 5th day, REL of E1 was significantly higher than that of M79 (P < 0.01), while the REL of M79 on the 6th day was the lowest among the three genotypes (20.87%) (Fig. 1C). The CAT activity of M79 was the highest (10.97) on the 6th day (Fig. 1D), which was significantly higher than that of E1 and H1. Although the WUE of all genotypes decreased after drought stress, M79 exhibited a significantly higher WUE on the 6th and 7th days (4.37 and 4.10, respectively) compared to E1 and H1 (Fig. 1E). These results demonstrated that M79 had better tolerance to drought stress than its parental lines as shown by phenotypical and physiological indexes.

RNA-seq data export, quality control, and sequence alignment

Leaves of M79 and its parental lines were sampled for transcriptomic sequencing six days after drought stress. A total of 18 libraries were constructed and sequenced using the HiSeq X Ten sequencing platform, and generated 8.55 × 10⁸ raw reads. After removing the linker and low-quality data, we obtained 8.17 × 10⁸ (95.57%) clean reads, consisting of about 122.58 Gb of clean data, and representing an average of 4.54 × 10⁷ clean reads, i.e., about 6.81 Gb of valid data per sample. Phred mass fraction Q30 (error rate 0.1%) ranged from 86.38 to 88.22%, with an average GC content of 56.30%. We aligned 93.13 to 94.32% of the valid data to the reference genome (Table S1). The FPKM density distribution (Fig. S1A) and FPKM box diagram (Fig. S1B) suggested that the density of the detected genes followed a standard normal distribution. These results indicated high quality and reasonable reproducibility of our sequence data.

Validation by qRT-PCR

To verify the reliability of our transcriptomic sequence data, we selected ten genes from all three lines for qRT-PCR, including genes encoding POD, No Apical Meristem, ATAF1/2, Cup-Shaped Cotyledon 2 (NAC) transcription factor, wax, lipid transfer proteins (LTPL78), Domain of unknown function (DUF538), expansin precursor, PsbP, Psb28, and two oxygen evolving enhancer proteins 3 (Table S2). Positive correlation coefficients between RNA-seq and qRT-PCR results was high and significant (R² = 0.975, P < 0.01, Fig. 2), indicating that the transcriptome sequencing results were accurate and reliable.

Comparative analysis of DEGs between M79 and its parental lines before and after drought stress

Comparative analysis of DEGs in non-stressed plants showed that M79 had 1359 and 648 genes up- and down-regulated, respectively, when compared to E1, and had 1496 and 1033, respectively, when compared to H1 (Fig. S2). To explore how these DEGs enhanced drought-resistance, a GO analysis was performed on highly DEGs identified between M79_CK and E1_CK, and between M79_CK and H1_CK. The DEGs found between M79 and E1 were mainly involved in ADP binding, RNA synthesis, post-translational modification, cell recognition, and carbohydrate metabolism (Table S3). In the comparison between M79 and H1, the DEGs were mainly related to protein kinase, iron binding, redox balance, and post-translational modification (Table S4).

We analyzed the DEGs between M79_DR and E1_DR, and between M79_DR and H1_DR. 5,258 (2739 up-regulated, 2519 down-regulated) DEGs were identified between M79_DR and E1_DR, and 3594 (1795 up-regulated, 1799 down-regulated) DEGs were identified between M79_DR and H1_DR (Fig. S3). GO analysis showed that the DEGs between M79_DR and E1_DR, and between M79_DR and H1_DR were significantly enriched in protein kinase activity, ATP binding, iron ion binding, carbohydrate metabolism, redox balance, and post-translational modification (Tables S5 and S6). KEGG analysis on these DEGs found 18 pathways (FDR < 0.05) after the comparison between M79_DR and E1_DR, including metabolism of glutathione, phenylalanine, porphyrin, chlorophyll, arginine, and proline, and biosynthesis of phenylpropanoid, carotenoids, flavonoids, cuticle, suberin, and wax (Fig. 3A). Comparison between M79_DR and H1_DR identified 15 different pathways (FDR < 0.05) including metabolism of glutathione, porphyrin and chlorophyll, and biosynthesis of carotenoids, brassinosteroids, cuticle, suberin and wax, and plant hormone signaling (Fig. 3B).

DEGs analysis of M79, E1 and H1 under drought stress

Compared to untreated plants, 3066, 1895, and 2148 DEGs were identified after drought treatment in M79, E1 and H1, respectively, with 1404, 1116, and 1328 up-regulated genes and 1662, 779 and 820 down-regulated genes in corresponding genotypes. Among these DEGs, 288 (208 up-regulated and 80 down-regulated) genes were expressed in all three genotypes, accounting for 9.39%, 15.20%, and 13.41% of all DEGs in drought-treated M79, E1 and H1, respectively (Fig. 4). GO analysis showed that these genes were significantly enriched in carbohydrate metabolism, iron ion and heme binding, oxidoreductase, POD, protein kinase activity, and plasma membrane osmoregulation (Tables S7 and S8). Among them, genes known to be involved in drought tolerance included two POD precursors (Seita.5G174100 and Seita.8G015200), two late embryogenesis abundant proteins (LEAs) (Seita.1G015800 and Seita.5G021400), and two aquaporins (Seita.3G082100 and Seita.1G264900). Furthermore, genes involved in photosynthesis, such as one early light-induced protein (Seita.2G053800), one pheophorbide a oxygenase (PaO) (Seita.1G348100), and one senescence-inducible chloroplast stay-green protein 1 (SGR1) (Seita.2G285600) (Table S9) were also found. In addition, many transcription factors including basic leucine zipper (bZIP), NAC, v-myb avian c viral oncogene homolog (MYB) and early responsive to dehydration (ERD) family members, as well as genes related to calmodulin, protein kinase, and hormone (gibberellin (GA), ethylene (ETH)) biosynthesis and signaling (Table S9).

Further analysis of DEGs that were specifically expressed before and after drought stress detected 908, 475 and 745 up-regulated genes, and 1368, 404, and 497 down-regulated genes in M79, E1 and H1, respectively (Fig. 4). GO analysis showed that the drought-specific DEGs in M79 were significantly enriched in the GO-terms associated with PSII oxygen-evolving complex, carbohydrate metabolism, redox balance, and iron ion binding (Table S10). In contrast, the drought-specific DEGs were mainly enriched in GO-terms associated with iron ion binding, redox balance, and nucleic acid and transcription factor activity in E1 (Table S11), and with iron ion binding, protein kinase activity, and fatty acid synthesis in H1 (Table S12). KEGG analysis showed that the DEGs in drought-stressed M79 were mostly involved in pathways such as phenylalanine biosynthesis, plant hormone signaling, porphyrin, chlorophyll metabolism, cuticle and wax biosynthesis, and arginine and proline metabolism (Table S13). The pathways included phenylalanine biosynthesis, phenylalanine metabolism, plant hormone signaling, linoleic acid metabolism, and glycerophospholipid metabolism in E1 (Table S14), and phenylalanine biosynthesis, plant hormone signaling, and carotenoid biosynthesis in H1 (Table S15). These results suggested that those DEGs enabled the three genotypes to respond differently to drought stress.

Expression and regulation of drought stress-responsive genes in M79

GO enrichment and KEGG pathway analysis indicated that these stress-induced DEGs in M79 were widely involved in signal transduction, transcriptional regulation, hormone signaling, redox regulation, osmotic regulation, photosynthesis, and other biological processes (Table S16).

Among the signal transduction-related genes that were differentially expressed in M79 after drought stress, 24 genes encode receptor kinases, of which 14 encoded wall-associated kinase receptor-like protein kinase (WAK-RLK) (eight up-regulated and six down-regulated), 27 genes encode protein kinases, of which 10 encode calcium/calmodulin dependent protein kinases (CAMK) (five up-regulated and five down-regulated), and six are Ca²⁺-related genes, including one up-regulated gene that encodes an EF hand family protein, and five genes encoding sodium/calcium exchanger (NCX) protein (three up-regulated and two down-regulated) (Table S16).

Ninety-six transcription factors were differentially expressed in M79 in response to drought stress, including 20 NAC, 19 APETALA2 (AP2), 14 transcription factors containing highly conserved protein domain (WRKY), seven bZIP, three ethylene response factor (ERF), and five dehydration-responsive element-binding (DREB) family transcription factors. These transcription factors played essential roles in M79 when it was exposed to drought stress (Table S16).

Many genes involved in phytohormone signaling were also responsible for drought tolerance in M79. We identified 17 DEGs that are related to auxin/indole-3-acetic acid (Aux/IAA) regulation, including nine of the OsIAA family (one up-regulated and eight down-regulated) and four belonging to the OsSAUR family (one up-regulated and three down-regulated). In addition, we identified one up-regulated cytokinin (CTK) dehydrogenase-related gene, ten DEGs involved in the GA regulation pathway, including six gibberellin 20 oxidases (four up-regulated and two down-regulated), three gibberellin 2-beta-dioxygenases (two up-regulated and one down-regulated), and one up-regulated gibberellin receptor GID1L2, and three DEGs related to ETH (two up-regulated and one down-regulated) (Table S16).

A total of 63 DEGs in M79 were identified to be involved in redox regulation, including genes encoding superoxide dismutase (SOD), glutathione peroxidase (GPx), oxidoreductase, POD, ascorbate peroxidase (APX), and lipoxygenase (LOX) (Table S16).

Forty-four DEGs related to osmotic regulation in M79. Among them, twelve DEGs were involved in proline metabolism, with ten up-regulated genes; seven DEGs were aquaporin genes with three up-regulated genes, and 12 DEGs belonging to ATP-binding cassette, subfamily G (ABCG) transporter family genes with six up-regulated genes (Table S16).

Forty-nine photosynthesis-related genes were differentially expressed in drought-treated M79 plants. These were involved in several photosynthetic metabolic pathways, including synthesis and degradation of chlorophyll, light energy absorption and transmission, PSII reaction center electron transfer, PSII reaction center electron transfer, and water oxidation. Among them, one PaO gene, one PsbP gene, one phytoene synthase gene, one scaffold protein in nitrogen fixation system (NifU) gene, and one ferrochelatase-2 gene were significantly up-regulated in M79 (Table S16). These genes maintained photosynthesis in M79 under drought.

Co-regulation analysis of drought-responsive DEGs in M79

The co-regulation study generated a regulatory network of 72 genes (Fig. 5), and the related genes were further divided into five groups. Genes in group A were mainly involved in signal transduction, including those encoding receptor kinases, Ca²⁺-related proteins, and protein kinases. Group B contained genes for GA, Aux/IAA, ETH, and CTK signaling. Group C consisted of transcriptional regulatory genes, including 4 bZIP, 1 DREB, and three WRKY family transcription factors. Genes in group D were drought-related, acting downstream of the molecular pathway responsible for drought tolerance in M79, including redox balance regulation genes (POD, GPx, and APX) and osmotic regulation-related genes (ion transporter, aquaporins, and proline synthesis-related genes). Group E contained photosynthesis-related genes, including two genes encoding oxygen evolving enhancer protein, three encoding PsbP proteins, and one encoding a ferrochelatase-2. The regulatory network was involved at all stages, and the drought response pathways in M79 may be essential for higher drought tolerance in M79 than its parental lines.

Responses of drought-treated M79 in photosynthesis-related pathways

As a C4 crop, the ability for millet to maintain photosynthesis under drought stress is an important indicator of drought resistance (Feng et al., 2015). The net photosynthetic rate of the three cultivars eventually decreased with prolonged drought stress. The net photosynthetic rate of E1 was the highest before drought treatment (16.82) in comparison with H1 (14.57) and M79 (13.77). After the 6th and 7th days of stresses, the net photosynthetic rate of M79 was the highest (11.52 and 10.67, respectively), reflecting the smallest decrease among the three cultivars in response to drought stress (Fig. 6A). Under drought stress, both Fv/Fm and ΦPSII showed a decreasing trend, and the parental lines E1 and H1 declined to a greater extent than M79 (Figs. 6B, 6C). These results demonstrated that the function of PSII was inhibited under drought stress, and that M79 maintained a relatively higher light energy utilization ratio than E1 and H1.

Functional annotation, GO enrichment and KEGG analysis identified 49 DEGs involved in the photosynthesis pathway responding to drought stress and showed higher expression levels in M79 than E1 and H1 (Fig. 6D). Further analysis of the photosynthetic pathway in M79 under drought stress showed that a gene encoding ferrochelatase-2 (Seita.4G016600) was up-regulated and involved in absorption and utilization of light energy during photosynthesis. Moreover, two genes encoding oxygen evolving enhancer protein (Seita.J002400 and Seita.1G208500) along with three PsbP genes (Seita.3G333700, Seita.9G561800 and Seita.5G442500) and one gene encoding the PSII reaction center Psb28 protein (Seita.5G446500) were involved in the breakage of water and oxygen release in the PSII reaction center. In addition, the NifU gene (Seita.8G058400) was up-regulated and involved in Fe-S cluster assembly in the PSII reaction center. Finally, two APX genes (Seita.7G023900 and Seita.9G444200) played a role in scavenging ROS (Fig. 6E). In addition, five genes (Seita.2G303000, Seita.5G234700, Seita.5G235200, Seita.6G024800 and Seita.9G545500) could be directly linked to the net photosynthetic rate and Fv/Fm based on the correlation analysis on physiological data and the expression of photosynthesis-related DEGs (Table S17). All of these genes regulated the photosynthetic pathway in M79 in response to drought stress.

The drought responsive pathway and related genes in millet

Previous studies on drought-responsive pathways indicate that drought-inducible genes from different varieties of the same crop likely play a relatively conserved role in their regulatory networks (Yamaguchi-Shinozaki & Shinozaki, 2006). However, drought-induced genes codify not only proteins that directly protect the cell structure and related metabolic pathways, but also regulators with roles in stress signaling, and forming a set of elements responding to environmental stress (Caldana et al., 2011; Vermeirssen et al., 2014). Among them, WAK and CAMK are central in plant responses to abiotic stress (Zhang et al., 2005; Lv et al., 2014). Protein phosphatase 2C (PP2C) belongs to a group of phosphatases involved in ABA signaling, and is a negative regulator of ABA signaling (Zhang & Gan, 2012). Among the DEGs that we identified from three drought-treated cultivars, four genes encoding CAMK kinases were up-regulated, three genes encoding WAK receptor kinases were down-regulated, and seven PP2C genes were up-regulated in all three cultivars (Table S9). Aux/IAA, ETH, and GA are phytohormones playing an important part in maintaining normal plant growth and development, and in reacting to abiotic stresses (Arraes et al., 2015; Gaion et al., 2017; Yu et al., 2017). Our study detected one down-regulated gene encoding for Aux/IAA, two ETH-encoding genes (one up-regulated and one down-regulated) and three GA-encoding genes (one up-regulated and two down-regulated) in all three drought-treated cultivars (Table S9). Transcription factor families such as NAC, MYB, bZIP, and basic helix-loop-helix (bHLH) also have a role in plant responses to both biotic and abiotic stresses (Katiyar et al., 2012; Puranik et al., 2013; Yong et al., 2014). Twelve DEGs belonging to these families were expressed in all three cultivars. They include five up-regulated genes encoding NACs, five genes encoding MYBs (three up-regulated and two down-regulated), an up-regulated gene encoding a bZIP and an up-regulated gene encoding a bHLH transcription factor (Table S9). Many protective proteins are also considered essential for protecting plants from damage caused by drought stress. For example, drought stress induces the expression of LEAs, which reduce water loss in plant cells and increase their WUE under adverse conditions (Wang et al., 2014). Genes encoding for protective proteins that were up-regulated in all three cultivars in our study included three coding for universal stress proteins (USP) and two for LEAs (Table S9).

A variety of intrinsic membrane proteins protect plant cells from abiotic stress by regulating the permeability of the plasma membrane (Kasim et al., 2015). Here, we detected two genes encoding aquaporins, one up-regulated and one down-regulated (Table S9). Plant cuticle is a hydrophobic protective layer that prevents water loss and protects plants from abiotic stress, such as those created from exposition to high temperature, drought, and salt (Ma et al., 2015). We detected both up-regulated genes encoding waxs in our experiments (Table S9). The antioxidant defense system in plants under drought stress is composed of ROS scavenging enzymes. Among them, CAT, SOD, APX, and GPx are essential to remove ROS and act synergistically to counteract oxidative damage caused by drought stress (Adriano et al., 2015). Our results indicated that one gene encoding POD was up-regulated and one down-regulated (Table S9). In addition, LTPs play an important role in response to biotic and abiotic stresses. We detected seven LTPL-coding genes up-regulated in all three cultivars, suggesting a positive role for these proteins in the drought defensive pathway (Safi et al., 2015 Table S9). All DEGs identified in drought-treated M79, E1 and H1 were involved in drought-defensive processes, and the transcription of many of them was up-regulated, suggesting that these genes play a positive regulatory role in drought response.

Molecular basis for better drought tolerance in M79 than its parental lines

RNA-seq analysis identified 5258 DEGs between M79 and E1, and 3594 between M79 and H1, indicating that the drought-tolerant cultivar M79 and its parental lines had different transcriptional profiles (Fig. S3). GO analysis of these DEGs showed that they were highly enriched in GO-terms such as membrane, protein kinase activity, transferase activity, carbohydrate metabolism, iron ion binding, ATP binding, heme binding, oxidoreductase activity, phosphorylation, and protein modification (Tables S5, S6). Among them, metal ion binding, electron-carrier activity, and expression of genes related to oxidoreductase synthesis increase the ability of plants to resist drought and high temperature (Rizhsky et al., 2004). In addition, phosphorylation is involved in stress responses in plants. In Arabidopsis, the SUCROSE NONFERMENTING1 (SNF1) kinase homologs 10 and 11 play an essential role in stress responses (Chen & Hoehenwarter, 2015).

KEGG analysis found 18 and 15 significant pathways that differentiate between M79_DR and E1_DR (FDR < 0.05), and between M79_DR and H1_DR (FDR < 0.05). These pathways include biosynthesis of secondary metabolites, plant-pathogen interaction, metabolic pathways, carotenoid biosynthesis, glutathione metabolism, amino sugar and nucleotide sugar metabolism, phenylalanine metabolism, phenylpropanoid biosynthesis, porphyrin and chlorophyll metabolism, diterpenoid biosynthesis, monoterpenoid biosynthesis, arginine and proline metabolism, glycerophospholipid metabolism, ubiquinone and other terpenoid-quinone biosynthesis, limonene and pinene degradation, flavonoid biosynthesis, cutin, suberine and wax biosynthesis, galactose metabolism, brassinosteroid biosynthesis, glycerolipid metabolism, plant hormone signal transduction, and alanine, aspartate and glutamate metabolism (Fig. 3). Among them, glutathione metabolism can reduce and eliminate oxidative damage caused by ROS, and it plays an important role in maintaining redox balance (Hicks et al., 2007). As an important osmotic regulator in plants, proline helps to maintain osmotic pressure, and it stabilizes proteins and cellular structures under drought stress (Vendruscolo et al., 2007). Phenylalanine and flavonoids also play an important role in adapting plants to stress and overcoming stress damage (Hernández, Alegre & Munné-Bosch, 2006; Babst et al., 2014; Pan et al., 2018). The cuticle and wax contribute in protecting plants from biotic and abiotic stresses, and in maintaining plant morphology (Pollard et al., 2008). Using high-throughput Illumina RNA-seq, Wang et al. (2016a) identified and quantified DEGs related to flavonoid biosynthesis in drought-stressed plants. Biosynthesis of brassinosteroids, monoterpenoids, porphyrins, chlorophyll, ubiquinone, and other terpenoid-quinones also is instrumental in adapting plants to drought stress and overcoming stress damage (Ksouri et al., 2016; Hajrah et al., 2017; Miao et al., 2017). In addition, KEGG pathway analysis indicated a significant enrichment in amino sugar and nucleotide sugar metabolism, and limonene and pinene degradation in response to abiotic stress (Singh et al., 2017; Zhang et al., 2017).

When challenged with drought stress, the transcriptional profile and ability to tolerate drought stress in M79 were significantly different from those of both parental lines. Crossbreeding is a means to promote recombination of parental genes, and is also a prerequisite for breeding superior offspring, that not only inherit desired traits from both parental lines, but also can exhibit novel F₁ phenotypes, heterosis and modified gene expression (Shivaprasad et al., 2012; Bell et al., 2013). However, the role of heterosis in defense mechanism of abiotic stress is still poorly understood (Dong et al., 2006; Korn et al., 2008; Singh, Sharma & Singh, 2010; Miller et al., 2015). Korn et al. (2008) showed that there is a strikingly strong correlation between heterosis, freezing tolerance, and flavonol content. The heterotic vigor for SOD, POD and CAT suggests an improvement of stress tolerance level in hybrids compared to the parental lines (Singh, Sharma & Singh, 2010). Miller et al. (2015) found that the levels of stress-responsive gene expression in parental lines could be used to predict biomass heterosis in hybrids. We therefore hypothesized that the reason for the higher drought stress tolerance in M79 than its parental lines was due to heterosis. However, this hypothesis needs further examination.

Molecular co-regulatory network for drought tolerance in millet

The millet genome has a complex molecular regulatory network to cope with drought, which can activate specific cell signaling pathways and induce transcriptional regulation that leads to enhanced cellular responses, increased expression of antioxidant-related genes, and accumulation of soluble substances (Lata, Sahu & Prasad, 2010). We analyzed DEGs in M79 before and after drought stress. Genes that are widely involved in biological processes such as signal transduction, hormonal signaling, transcriptional regulation, redox regulation, osmotic adjustment, and photosynthesis formed a drought-tolerance regulatory network (Table S16).

Many signal transduction-related genes were up-regulated after drought treatment, including five genes encoding NCX proteins (Table S16). This gene family is involved in Ca²⁺ signaling, and dehydration induces the expression of OsNCX3, OsNCX10, and OsNCX15 in rice (Singh et al., 2015). Therefore, NCXs in M79 may be involved in signal transduction during drought stress to activate the expression of downstream genes. In addition, we found 14 DEGs (eight up-regulated) encoding WAKs, and 10 (five up-regulated) encoding CAMKs (Table S16). WAKs belong to a receptor-like kinase gene family, and respond to abiotic stress acting on the signal transduction between the cell wall and cytoplasm (Zhang et al., 2005), whereas CAMKs have roles in Ca²⁺ signaling and protein phosphorylation (Chen, Zhang & Cheng, 2017). These up-regulated genes positively regulated drought-related signaling pathways in M79, contributing to drought stress tolerance with a rapid activation of downstream gene expression.

ETH, CTK, Aux/IAA, and GA are all involved in plant response to abiotic stresses (Werner & Schmülling, 2009; De Diego et al., 2013; Colebrook et al., 2014; Zwack & Rashotte, 2015; Yang et al., 2015; Yu et al., 2017). CTK and Aux/IAA negatively regulate ABA-induced stomatal closure (Werner & Schmülling, 2009; De Diego et al., 2013). Colebrook et al. (2014) showed that decreased GA content and transcriptional alteration of related genes inhibit plant growth and development under various abiotic stress conditions. Yang et al. (2015) demonstrated that ETH is involved in salt stress-related responses in rice, and that it plays a vital role in regulating biotic and abiotic stress responses. In this study, we found two, one, three, and seven up-regulated genes involved in ETH, CTK, Aux/IAA, and GA pathways, respectively, in M79 following drought stress (Table S16). These genes maintain growth and development in plants under drought stress by regulating hormone balance.

Transcription factors play a central role in biotic and abiotic stress responses, and in the regulation of various biological processes. AP2/DREB, WRKYs, ERF, bHLHs, bZIP and NAC transcription factor families are involved in transcriptional regulation in response to stress (Lata et al., 2014; Li et al., 2014; Muthamilarasan et al., 2015). Expression analysis of the millet AP2/ERF genes SiAP2/ERF-069, SiAP2/ERF-103, and SiAP2/ERF-120 showed that they were all up-regulated under drought stress and therefore may play a positive role in this process (Lata et al., 2014). SiARDP belongs to the DREB family of transcription factors and is one of the target genes of SiAREB; it participates in the ABA-dependent signaling pathway. Overexpression of SiARDP improves drought resistance in millet (Li et al., 2014). In addition, Muthamilarasan et al. (2015) showed that the SiWRKY genes SiWRKY066 and SiWRKY08 give an essential contribution in response to abiotic stress. In this study, multiple members of these families were differentially expressed after drought stress, including 20 NACs, 19 AP2s, 14 WRKYs, seven bZIPs, three ERFs, and five DREBs, indicating their involvement in the responses to drought stress (Table S16). In addition, Komivi, Diaga & Ndiaga (2016) showed that 90% of heat shock factors (HSFs) respond to drought stress in sesame seedlings, and that two HSF transcription factors are significantly up-regulated after drought stress, suggesting that these genes might contribute to this process (Table S16).

Upon drought stress, gene expression is pivotal in protecting plants from oxidative damage. For example, SOD converts superoxide into the less toxic H₂O₂, which is then reduced to H₂O by POD, APX, and GPx (Alscher, Erturk & Heath, 2002; Fecht-Christoffers et al., 2006; Islam, Manna & Reddy, 2015). In our study, drought stress up-regulated 10 genes encoding PODs, two genes encoding GPx, and one gene encoding APX in M79. These genes function as positive regulators removing ROS and maintaining a redox balance (Table S16).

Our study detected a large number of genes related to osmoregulation in M79 after drought stress, including 12 (10 up-regulated) involved in proline metabolism, 12 (six up-regulated) ABCG transporters, seven (three up-regulated) aquaporins, and two (both up-regulated) ATP-binding cassette (ABC) transporters (Table S16). Proline is an important macromolecule serving as an osmotic regulator in stress-defensive responses, since its accumulation can relieve damage caused by osmotic stress under drought (An et al., 2013). ABC transporters utilize ATP hydrolysis to transport osmotic regulators such as amino acids, peptides, carbohydrates, lipids, hormones and metal ions (Jeong et al., 2014). Aquaporin proteins assist the plant to combat abiotic stress by regulating the permeability of the plasma membrane (Kasim et al., 2015). ABCG transporters help to biosynthesize protective cuticles and wax, transporting lipids or to regulate phytohormone homeostasis transporting indole butyric acid and ABA (Yadav et al., 2014). Expression of the above genes can regulate the osmotic potential of M79 cells under drought stress to reduce injury.

Co-regulation analysis of these DEGs in M79 revealed a regulatory network consisting of 72 genes, which might contribute to the excellent drought resistance of M79. This system includes signal perception and transduction, hormone signaling pathways, transcriptional regulatory factors, and downstream functional genes (including ROS removal factors, ion transporters and osmotic regulators) (Fig. 5). Although the results of our co-regulation analysis remain to be further verified, the co-regulation network provides an important theoretical basis to propose a model for the molecular mechanisms of drought tolerance in millet.

Maintenance of a high photosynthetic rate is an important indicator of drought tolerance in crop plants

Photosynthesis is the basic metabolism regulating crop growth and final yield. The maintenance of photosynthetic rates under drought stress is essential for drought tolerance in crops (Galmés, Medrano & Flexas, 2007; Chaves, Flexas & Pinheiro, 2009). We found that photosynthetic rate and light energy utilization in E1 and H1 was significantly lower than in M79 after drought stress, suggesting that M79 can maintain higher photosynthesis under drought (Figs. 6A– 6C). Therefore, photosynthetic rate under drought stress is not only related to photosynthetic capacity, but also to drought tolerance (Zhang, Li & Xiao, 2016).

Photosynthesis in chloroplasts converts light into chemical energy that is used for plant growth and development. Under drought stress, O₂ produced in chloroplasts can receive electrons from the photosynthetic electron transport chain to become O⁻², which can cause oxidative damage to photosynthetic pigments and the plasma membranes (Gill & Tuteja, 2010; Exposito-Rodriguez et al., 2017). Previous studies showed that osmoregulation and antioxidant capacity of plants contribute to maintaining photosynthetic capacity (Ramachandra Reddy, Chaitanya & Vivekanandan, 2004). Our co-regulatory analysis showed that osmotic regulation and antioxidation played a vital role in the management of drought tolerance in M79 (Fig. 5; Table S16), explaining why M79 was able to maintain a high net photosynthetic rate under drought stress.

Our results identified 49 photosynthesis-related DEGs in drought-treated M79. Among them, 11 genes were involved in the photosynthetic pathway, including one gene encoding ferrochelatase-2 (Seita.4G016600), two genes encoding oxygen evolving enhancer proteins (Seita.J002400 and Seita.1G208500), three PsbP-encoding genes (Seita.3G333700, Seita.9G561800 and Seita.5G442500), one gene encoding PSII reaction center Psb28 protein (Seita.5G446500), one NifU-encoding gene (Seita.8G058400), and two APX-encoding genes (Seita.7G023900 and Seita.9G444200) (Fig. 6E). Previous studies showed that ferrochelatase is related to the absorption of light by the light-harvesting complex (LHC) proteins (Suzuki et al., 2002; Espinas et al., 2016). PsbP and Psb28, two subunits of the PSII reaction center, are involved in the water photolysis and oxygen release during photosynthesis (Mabbitt, Wilbanks & Eaton-Rye, 2014). NifU plays an important role in the synthesis and assembly of the Fe-S cluster in PSI (Yabe et al., 2004). Under drought stress, electron transport and photosynthetic phosphorylation in chloroplasts produce a large amount of ROS, while APX effectively removes them and reduces oxidative damage to plants (Jiang et al., 2016; Exposito-Rodriguez et al., 2017). Up-regulation of these genes enables M79 to maintain a relatively high level of photosynthesis under drought stress, and to resist the damage caused by drought. Through the analysis of the photosynthetic pathway in drought-stressed M79, we provide an initial insight to the molecular mechanism of drought resistance.