JIA-2019-11

2461 XU Bing-qin et al. Journal of Integrative Agriculture 2019, 18(11): 2457–2471 0 h control, then it increased gradually by 40.00–42.17% in the leaves from 48 to 72 h of drought treatment (Fig. 1-E and Appendix B). In addition, JA content increased greatly in the leaves of the two cultivars under 24 h drought treatment, with an increase of 621.72% in DM and 198.81% in HN after 24 h drought stress (Fig. 1-F andAppendix B). Then after 48 and 72 h drought stresses, the JAof two cultivars suffered a decline to various degrees compared with the level at 24 h. The dynamic physiochemical properties indicated that the time-course of drought stress caused rapid, but differential degrees of damage in leaves of the two foxtail millet cultivars. 3.2. RNA-seq analysis of two cultivars of foxtail millet under drought stress Based on the physiological changes that occurred, especially MDA, ABA, and JA contents of the two cultivars above under drought treatments, the 24 h drought treatment is considered as an important point at which the plants had sensed the drought and responded quickly. So RNA isolated from plants of the two cultivars that had been subjected to drought stress for 0 and 24 h, with each condition having three biological replicates, were selected to construct 12 cDNA libraries denoted as: DM1 for DM in 0 h; DM2 for DM in 24 h; HN1 for HN in 0 h; and HN2 for HN in 24 h. As a result, 0.35 billion sequence reads were generated in total, with 150 bp for each read. After trimming adapters, eliminating low-quality reads, and removing reads that contain rRNA with Bowtie version 2.2.4, 0.30 billion clean reads (98.65% of the generated data) with approximately 28 million reads from each sample were recorded, then mapped to the reference genome of foxtail millet using TopHat2 (Kim et al . 2013). As a result, among the high-quality clean reads generated from the 12 samples, uniquely mapped reads were 82.99 to 85.70%, while total mapped reads were 85.92 to 86.99% (Table 1). 3.3. Differential expression analysis of transcripts Across the four samples, an average of 26 406 (73.7%) known protein-coding genes and 944 new genes were detected (Appendix C). Data quality assessments and correlation analyses between samples were conducted using heatmap visualization of Pearson’s r (Stigler 1989) of expression gene amounts in any two samples and principle component analysis (PCA). Hierarchical clustering analysis of all samples was conducted using R Project version 3.3.3 (https://www.r-project.org/) (Appendix D), showing that the repeatability of the three biological replicates was satisfactory. The heatmap based on the Pearson’s r (Appendix D) suggested that drought stress had greater effects than genotype on the gene expression patterns, and that interactions existed between drought treatment and genotypes. This agreed with the report of Tang et al . (2017), but the interaction in this report had less influence. Differentially expressed genes (DEGs) between the two groups, DM1/DM2 and HN1/HN2, were found with edgeR version 3 (http://www.bioconductor.org/packages/release/ bioc/html/edgeR.html) with an FDR of 0.05, |log 2 FC|>1. A total of 5 554 DEGs detected were either up- or down- regulated (Fig. 2-A). Of those DEGs for DM, 1 644 were up-regulated under 24 h drought stress (DM2) and 1 870 were down-regulated compared with the 0 h (DM1); while for HN, 1 311 DEGs were up-regulated under 24 h water- deficit stress (HN2), and 729 genes were down-regulated compared with the 0 h (HN1) (Fig. 2-A). Volcano plot was used to visually display the expression difference and screening of DEGs using R Project (https://www.r-project. org/) (Fig. 2-B). 3.4. GO and KEGG analysis of differentially expressed genes (DEGs) Based on GO annotation, the drought stress mainly affected three biological processes: single-organism processes, metabolic processes, and cellular processes. For molecular functions, drought stress affected binding and catalytic activity, and for the cellular components, it affected cells, cell parts, and organelles (Fig. 3). The main changes of molecular functions and cell components in HN1<HN2(+) were similar to those of DM1<DM2(+). Based on the KEGG database, these DEGs were significantly enriched in metabolic pathways, especially amino acid metabolism and carbohydrate metabolism in DM1<DM2(+) and HN1<HN2(+) (Fig. 4), by performing pathway enrichment analysis of the DEGs. Plant hormone signal transduction also made up a large proportion, with 47 genes in DM1<DM2(+) and 36 genes in HN1<HN2(+) (Fig. 4). Interestingly, the KEGG terms “starch and sucrose metabolism”, “galactose metabolism”, and “glycolysis/ gluconeogenesis” were only enriched in HN with respect to drought-induced DEGs, which might contribute to the high soluble sugar content in HN after 24 h drought stress (Fig. 4). Genes that participated in “ascorbate and aldarate metabolism” and “glutathione metabolism” are particularly enriched in the genotype HN, and they may participate in the ascorbate-glutathione (AsA-GSH) cycle involved in ROS to reduce the MDA level (Fig. 1). 3.5. Cross-comparison of differentially expressed genes in the two pairs DM1/DM2 and HN1/HN2 In order to gain further insight into the regulatory network