Effect of drought stress on physiological indicators in Axonopus compressus
In Figure 1, it was revealed that the rate of H2O2 (Figure 1A) production and MDA concentration (Figure 1B) increased in drought treated plants. The MDA values and H2O2 production rate increase with the prolonged drought stress treatment, whereas MDA values and the production rate of H2O2 at 36 h after drought induction were significantly increased (p < 0.01). Similarly, the leaf water potential and relative water content decrease significantly under drought, contrarily, the electrolyte leakage increased over the drought period in Axonopus compressus. As in the case of MDA and H2O2, the values of EL were also increased as the time of drought stress increase. Whilst, the values of leaf water potential and relative water content decrease with the prolonged drought stress (Figure 1C).
RNA-Seq of Axonopus compressusand de novo assembly
RNA-Seq of the six libraries of cDNA (DS and CK with three replicates each) resulted in 352.84 million raw reads, of which approximately 347.73 million clean reads de-novo compiled into contigs utilizing Trinity software, tend to range from 54.90 million to 62.25 million reads for each library and approximately 97% of reads with a quality score of Q20 (99 % accuracy) (Additional file 1; Table S2). The contigs were assembled into 278042 unigenes with an average length of 1066 bp and an N50 length of 1430 bp. All unigenes were over 200 bp long, of which 11.77 % (32718 unigenes) were 2,000 bp long. We evaluated the assembled data of transcriptome by employing Benchmarking Universal Single-Copy Orthologs (BUSCO) package. The number of total BUSCOs searched was 1440. The number of complete BUSCOs in the combined transcriptome was 1080 (75.00%), and only 198 core genes (13.75%) are fragmented in Axonopus compressus (Additional file 2; Figure S1).
Functional annotation and classification of unigenes
The assembly of unigenes from Axonopus compressus was annotated by BLASTX (E-value <10−5) in different databases. A total of 135716 (48.81%), 188334 (67.64%), 120455 (43.32%), 189422 (68.13%), 100588 (36.18%), 153958 (55.37%) and 128626 (46.26%) unigenes would have significant levels (E-value <10−5) against KEGG, NR, Swissprot, Trembl, KOG, GO, and Pfam, respectively (Figure 2A). Out of 278042 unique of high-quality sequences, 191893 (69.03%) unigenes complemented a sequence significantly in at least one of the seven public databases, The five key online databases (GO, KEGG, KOG, NR, and Trembl) were picked out of seven databases to draw a Venn diagram (Figure 2B), the unigenes with significant values (E-value10-5) are also noticed at each junction of the Venn diagram, in which 94 unigenes corresponds to five databases.
The alignment of sequence homology revealed that 39439 (20.94%) sequences that were found in resemblance with Setaria italica; 35283 (18.73%) sequences had significant hits for Zea mays, followed by Sorghum bicolor (35034, 18.6%), Penicum hallii (23894, 17.47%), Dichanthelium oligosanthes (18896, 10.03%) and Oryza sativa japonica group (6551, 3.48%). Brachypodium dischyon which included 1675 (0.89) and 9.86% of the sequences (18562) were homologous to other species (Figure 2C). The expression of genes has biological variability between different individuals, and there are differences in the degree of expression between genes. The reads per kilobase of exon model per million of aligned reads (FPKM) values were calculated as normalized expression estimates for each gene model in each sample. To evaluate the major spectral variance between CK and DS treatment samples, we performed PCA (Figure 2D). Results showed a distinct separation of CK to DS sample treatments, which reveals that both groups of the treatment samples exhibit different spectral positions on the PCA chart (Figure 2D).
To demonstrate the accuracy of the unigenes prediction, the nucleotide sequences from the transcriptome of Axonopus compressus were used as a query in a BLAST search with a threshold E-value <10-5. The results showed 3,300 sequences of Axonopus compressus in the NCBI database. This analysis renders an assessment for the consistency of unigenes sequences in the current dataset. Based on the high-score BLASTx matches in the GO proteins database, The BLASTx high-score Predicated matches in the computer database of Go proteins were confirmed, and a total of 153958 unigenes were categorized with Blast2GO (E-value <10-5) and were designated at least once in GO. As shown in Fig 3, the unigenes referred to three main categories of GO and 59 subcategories, namely biological processes (BP), with 28 main sub-classes (405025 unigenes); cellular compartments (CC), with 18 main sub-classes (460094 unigenes); and molecular functions (MF), with 13 main sub-classes (194989 unigenes).
Cellular processes (21.43%) were the biggest subgroups in the category of biological processes, metabolic process (19.70%), biological regulation (9.24%), response to stimulus (9.22%), and regulation of biological process (8.36%). The largest subgroups in the cellular component category were cell (22.53 %) followed by cell part (22.47 %), organelle (17.22%), and membranes (11.72%) respectively. Similarly, the main subgroups in the molecular function category were binding and catalytic activity, which contribute 77.78% and 69.28% respectively, and 115199 unigenes associated with molecular function. Within the Axonopus compressus unigenes, 112492 (42.64%) were classified (E-value <10-5) into twenty-six KOG clusters (Figure 3B). The biggest groups which includes; 1) general function prediction, only (21955 genes, 19.52%); 2) posttranslational modification, protein turnover, and chaperones (11636 genes, 10.34%); 3) signal transduction mechanisms (11302 genes, 10.05%); 4) Function unknown (5785 genes, 5.14%); and 5) carbohydrate transport and metabolism (5616 genes, 4.99%).
Metabolic pathway analysis of Axonopus compressus by KEGG
A total of 135717 (51.44%) of the 263835 Axonopus compressus unigenes possessed significant correspondence in KO. All these unigenes have been restricted to 103 KEGG pathways in 5 major categories (Figure 3B). The pathways of KEGG (which includes 4303 unigenes) were the members of a major group, metabolism (D), 921 associated to genetic information processing (C), 184 related to cellular processes (A), 347 involved to environmental information (B), and 151 related to category (E) of organism systems (Additional file 2; Figure S2).
CDS prediction in Axonopus compressus
The BLASTx protein database (NR and SwissProt database) identified 278042 unigenous CDSs, of which 23652 unigenes were larger than 500 bp, 12520 unigenes were larger than 1,000 bp, and 32718 unigenes were larger than 2,000 bp. In addition, 43,713 unigenes were not linked to the NR and SwissProt database systems. The Estscan (Version; 3.0.3) software was used to interpret their ORF, frequency distribution, length, and related amino acid sequences of the unigene CDSs.
Differentially expressed genes (DEGs) analysisof Axonopus compressus
Amongst the differentially expressed unigenes, the expression of 7444 differs substantially between samples treated with drought-stress (DS) and control (CK) samples. Under drought treatment, 4249 numbers were up-regulated and 3195 numbers were down-regulated (p<0.05) (Figure 4A). The expression profiles of DEGs were also presented through cluster analysis that showed significantly different responses of CK and DS treatments in the Axonopus compressus (Figure 4B).
GO pathway enrichment analysis
By employing the Gene Ontology (GO) and the DEG enrichment analysis in Axonopus compressus, 23766 DEGs were categorized into three GO groups and 4592 numbers associated the DS with CK (number of DEGs annotated in more than one term), out of which 2766 items that were associated to BP, 477 were linked with CC and 1349 were related to MF.
The key 50 DEGs dramatically enriched in three GO categories are presented in Additional file 2; Figure S4. In the top 50 significantly enriched identities (Corrected P-Value<0.05), 18 were found to be linked with BP [carotenoid biosynthetic process GO:0016117, carotenoid metabolic process GO:0016116, cellular amine metabolic process GO:0044106, cellular biogenic amine metabolic process GO:0006576, cellular response to heat GO:0034605, cellular transition metal ion homeostasis GO:0046916, glutamine family amino acid biosynthetic process GO:0009084, Group II intron splicing GO:0000373, heat acclimation GO:0010286, oligosaccharide catabolic process GO:0009313, raffinose catabolic process GO:0034484, raffinose metabolic process GO:0033530, regulation of seed germination GO:0010029, regulation of seedling development GO:1900140, response to high light intensity GO:0009644, response to hydrogen peroxide GO:0042542, tetraterpenoid biosynthetic process GO:0016109, tetraterpenoid metabolic process GO:0016108], and 6 items were related to CC [chloroplast nucleoid GO:0042644, DNA packaging complex GO:0044815, Nucleoid GO:0009295, Nucleosome GO:0000786, plastid nucleoid GO:0042646, protein-DNA complex GO:0032993], and 14 items were related to MF [4-coumarate-CoA ligase activity GO:0016207, alpha-galactosidase activity GO:0004557, arogenate dehydratase activity GO:0047769, carbon-nitrogen ligase activity, with glutamine as amido-N-donor GO:0016884, ferric iron binding GO:0008199, ferroxidase activity GO:0004322, galactosidase activity GO:0015925, oxidoreductase activity, acting on single donors with incorporation of molecular oxygen GO:0016701, oxidoreductase response, acting on single donors with molecular oxygen incorporation, incorporation of two atoms of oxygen GO:0016702, oxidoreductase response, oxidizing metal ions, oxygen as acceptor GO:0016724, raffinose alpha-galactosidase activity GO:0052692, water channel activity GO:0015250, water transmembrane transporter activity GO:0005372] (Additional file 3).
Axonopus compressus KEGG pathway enrichment analysis
The KEGG pathway enrichment analysis for DEGs revealed that 2569 DEGs participated in 143 various types of pathways in Axonopus compressus. By comparing DS with CK, 2747 DEGs were found to be up-regulated and 2502 DEGs have been identified as down-regulated in deficit water. The twenty key pathways which were found significantly enriched by comparing DS with CK are displayed in (Figure 3C). The pathways includes Valine, leucine and isoleucine biosynthesis (215 unigenes), Stilbenoid, diarylheptanoid and gingerol biosynthesis (330 unigenes), Starch and sucrose metabolism (1642 unigenes), RNA degradation (1536 unigenes), Protein processing in endoplasmic reticulum (1628 unigenes), Porphyrin and chlorophyll metabolism (559 unigenes), Plant hormone signal transduction (1672 unigenes), Phenylalanine, tyrosine and tryptophan biosynthesis (601 unigenes), Metabolic pathways (1659 unigenes), MAPK signaling pathway plant (1661 unigenes), Linoleic acid metabolism (323 unigenes), Glyoxylate and dicarboxylate metabolism (1050 unigenes), Glycosphingolipid biosynthesis - globo and isoglobo series (103 unigenes), Galactose metabolism (735 unigenes), Carotenoid biosynthesis (506 unigenes), Carbon fixation in photosynthetic organisms (820 unigenes), Biosynthesis of secondary metabolites - unclassified (223 unigenes), Biosynthesis of secondary metabolites (1661 unigenes), Biosynthesis of amino acids (1655 unigenes), Alanine, aspartate and glutamate metabolism (639 unigenes) (Additional file 2; Figure S3).
Drought stress associated differentially expressed transcription factors (TFs)
The transcription factors play a crucial role in the growth and development of plants and therefore can stimulate and/or suppress transcriptional gene expression to sustain normal physiological functions in deficit water stress. In the current study, 352 TFs were observed in Axonopus compressus based on iTAK, including 270 transcription factors and 81 transcription regulatory factors that were associated with 32 and 16 families, respectively (Additional file 1; Table S3, S4). The differential analysis revealed that 270 transcription factors were different as shown in the comparison between DS and CK, 216 transcription factors have been noted to be up-regulated while 54 transcription factors were found down-regulated.
qRT-PCR Validation of Transcripts in Axonopus compressus under Drought
To check the validity of RNA-seq results, the expression of six transcripts for the three biological replicates was analyzed by qRT-PCR (Figure 5). All the genes were identified as significantly different expression values (DEGs). The findings of qRT-PCR revealed that the pattern of expression of these genes was identical with the findings in the RNA-seq analysis. To validate the results, the six drought-responsive genes (NAC, MAP Kinase1, MYB2, PIP1, WRKY1, ABI5) were verified using qRT-PCR. The different genes were amplified to compare the gene expression and RNA-seq results, and the findings of this study revealed that the expression of drought-responsive genes supporting the results of the RNA-seq.