Functional annotations of ESTs of Stevia rebaudiana involved in abiotic stress signaling through computational approach.

Stevia rebaudiana (Bertoni) is known as a natural sweetener plant to produced steviol glycosides. The steviol glycosides biosynthesis is limited in S. rebaudiana plants due to the alteration in the environmental circumstances such as drought, cold, salt and light. These environmental circumstances are a common side-effect in plants affecting the plant growth, metabolism and yield of secondary metabolites. Due to absence of complete genome annotations, the plant metabolites signaling is difficult in order to get the exact enzymatic flow to the product. In this article, we have analyzed the ESTs of S. rebaudiana and predicted their role in plant signaling in term of cis-regulatory elements, their biological function respect to abiotic stress. Further, the predicted abiotic stress responsive factors were also analyzed in order to predict the relevant genes or proteins function in comparison with the genome of Arabidopsis thaliana. Total 5,548 ESTs of stevia were retrieved from NCBI database. EST-contigs assembled from 5393 were 619 contigs and 2,894 singletons elements were identified by assembler program. Due to short expressed sequences related to singletons, it is excluded for further study. Further, retrieved ESTs were resulted in to 619 EST-contigs by using the clustering method. Out of 619, 15 contigs belongs to transcription factor families while 292 contigs, belongs to five enzyme classes. Out of 619 contigs, the 529 contigs showing the correct gene ontology in term of biological process (BP), molecular Functions (MF) and cellular component (CC). Further, these contigs were also screened for metabolic pathways analyses using KEG database. In this, 390 metabolic pathways and 67 involved for signal transduction were identified from 619 contigs. The co-expression analysis was revealed by gene investigators and STRING 10.0 with 0.40 correlations and 0.9 mutual connection. In projected PPI network, the recognized factors (WD40-like protein, MYB-HB like, AP2-EREBP, C2H2, Hap3/NF-YB, bHLH, C2C2-CO-like, CW-Zn, FHA-SMAD, Nin like, SBP3, TIFY, Tc-PD, Znf-B and bIP) belong to plant signaling and MAPK signaling pathways. These TFs introduce as a candidate genes responsive factors may lead to enhanced plant growth and metabolism by overexpression.

Introductions

Stevia rebaudiana (Bertoni; Asteraceae) is a medically important medicinal plant known as natural sweetener generally found in worldwide (Gerami et al., 2020). The plant has high sweetening compounds eg. steviol glycosides (diterpene compound) that are accumulated in the leaves (Lucho et al., 2019). Steviol glycosides (stevioside and rebaudioside) are the most commonly used in the control of diabetes across the globe (Gerami et al., 2020, Kumar et al., 2013, Richman et al., 2005).

The environmental stress, especially abiotic stress has serious concern which effect on plants health and reduced crop productivity (Fahad et al., 2017). These stresses include, salinity, drought, cold etc. and adversely affect the growth, development and plant metabolism (Wu et al., 2018, Acosta-Motos et al., 2016, Rao and Ravishankar, 2002, Hasegawa et al., 2000). The abiotic stress responses downregulate the plant growth, affect the cellular pools i.e., molecular and biochemical changes in plants (Sheshadri et al., 2016). In recent past, some researchers have described abiotic stress responsive gene (s) and cis-regulatory elements (TFs) especially those associated with transcriptional regulation of associated genes to play an important role in plant growth and metabolism (Sheshadri et al., 2016, Joshi et al., 2016). The abiotic stress responsive elements bind to the specific recognition sites of its genes and affects enzymatic pools in the metabolic process of plants (Duraisamy et al., 2016, Bussemaker et al., 2001). More recently, some of these stress responsive elements or factors such as MYB, bHLH, WRKY, NAC, ZIP9, ABF3 and ZIP1 are characterized and targeted as transcriptional regulators for plant development and metabolism under abiotic stress conditions (Chen et al., 2016, Tang et al., 2012, Ang et al., 1998).

To resolve this issue and to understand how transcriptional factors are playing an important role in biosynthesis of metabolites is still a herculean task. This due to the fact that molecular mechanisms and its multi-dimensional network i.e., genome organization is still completely unknown. The regulation of transcriptional factors is not a single step rather it involves post-transcriptional, translational and post-translational regulation also. To enhance the capabilities of plant growth and metabolism, researchers must target enzymatic pathways and regulations of the particular genes or TFs (transcription factors). However, the complete annotation of genome or ESTS though sequencing or using bioinformatics tools can help to gain more knowledge to understand the complexity of pathways and TFs role (Alam and Balawi, 2021, Yang et al., 2016). More recently, Singh et al. (2017), pointed out that some genes e.g., CYPs, UGTs, HMGR, DXS and KA13H genes, MYB, WRKY, bHLH, TFs played an important role in steviol glycosides production.

In this investigation, we have emphasized that the regulatory elements are the key factors that are vital in stevioside and rebaudioside A biosynthesis under the stress conditions. Stevia rebaudiana ESTs were retrieved from NCBI and role of cis-elements responsive genes was predicted in term of transcription elements, gene ontology (GO) enrichment analysis by comparing the genome of Arabidopsis thaliana.

Materials and methods

Data retrieval and assembling

5,548 ESTs of stevia retrieved from NCBI data base were masked to remove inappropriate sequence clustering. The ESTs were then masked for vector sequence, genomic repeats, sequencing artifacts, and low complexity sequence (including poly-A tails) using EGassembler (https://www.genome.jp/tools/egassembler/; Masoudi-Nejad et al., 2006). The retrieved sequences of ESTs obtained from NCBI database were assembled into clusters, constructed on the basis of sequences resemblance using CAP3online program (http://doua.prabi.fr/software/cap3,), with default (Huang and Madan, 1999). Each cluster were separately analyzed and identify the overlap and cut off range on the basis of minimum sum of nucleotides. The sequences that exhibited low similarity with other ESTs, could not be grouped and resulted in singletons.

Functional annotations

The functional characterization of assembled of EST-contigs was further analyzed by Blast2GO program, translated in all reading frames for the identification of potential translation products using BLASTX through Omicsbox full package software (Blast2GO program; www.biobam.com/omicsbox; Conesa and Götz, 2008). The assembled contigs were further used for the annotation of contigs against the proteins database for their functional characterization.

Abiotic stress signaling elements and Enzymes and KEG pathways

The assembled EST-contigs were further analyzed through PlantTFcat tool (http://plantgrn.noble.org/PlantTFcat; Dai et al., 2013) for TFs, similarly BLAST2GO program used for enzymes analyses (www.biobam.com/omicsbox). The SBH (single-directional best hit) process was applied for KEGG analysis to predict the enzymes pathways through Automatic Annotation Server (KAAS) 60, (http://www.genome.jp/kegg/kaas; Moriya et al., 2007).

Gene ontology network analyses

As per GO vocabulary, the EST-contigs sequences were characterized by OMICSBOX (www.biobam.com/omicsbox) to predoc the role of contigs in biological functions (BP, MF, and CC). Gene Ontology (GO) enrichment analysis for stress-responsive genes were carried out to determine their role in biological processes, molecular function, and cellular component through OMICS BOX (www.biobam.com/omicsbox; Conesa and Götz, 2008), and while gene ontology network or co-expression analyses of predicted TFs were analyzed by STRING 10.0 (http:// string-db.org; Szklarczyk et al., 2014) comparing Arabidopsis thaliana genome.

Results

EST-Contigs assembly

5,548 ESTs retrieved from NCBI db (database) of Stevia rebaudiana (Bertoni) were analyzed and the total length of 3162704 bp with 42.55% GC count were obtained (Table 1). The approx. length of predicted contigs from these ESTs was obtained as 620 bp. The entire essentials were screened and clipped into 17116 bp that constituted 0.54% of the entire size of EST-sequence, remaining EST sequences assembled into 619 contigs and 2,894 singletons using the same program in EGassembler. Out of 5,548 ESTs, only 11.47% accounted for assembled into ESTs contigs. The singletons showing partially expression signals could not assembled not further considered for analysis as they may be result from contamination.

Table 1

EGassembler based summary of Stevia rebaudiana ESTs.

Features	Numerals
EST sequences	5,548
Total number of EST nucleotides	3162704 bp
Clusters count	5,393
Singleton	2894
Contigs	619
Average GC content	42.55%

Identification of TFs and enzymes

A total of 23 contigs were screened from 619 contigs resulted belongs to the 15 putative TF families. These putative factors belong to WD40-like was found to be the most abundant category (21.43%), followed at MYB-HB like, AP2-EREBP, C2H2, Hap3/NF-YB, bHLH (14.29 each) while C2C2-CO-like, CW-Zn, FHA-SMAD, Nin like, SBP3, TIFY, Tc-PD, Znf-B and bIP (7.14% each; Table S1, Fig. 1). Similarly, 292 contigs were classified and characterized into 231 sub enzymes belongs to the six main classes of enzyme classes; hydrolases (35.49%) observed as abundant followed by oxidoreductases (32.46%), transferases (27.70) and lyases (17.31%) as compared to other while isomerases (6.91%), ligases (3.89%) and translocases (2.6%) (Fig. 2, Table S2).

Fig. 1

TFs distributions from the assembled EST-contigs of Stevia rebaudiana.

Fig. 2

Six main classes of enzymes belong to the assembled EST-contigs of Stevia rebaudiana.

Functional characterization of developed EST-contigs

The 619 predicted contigs were further analyzed by using the BLASTx program to translated the EST-contigs against the nucleotide database of NCBI. In this study, more than 50% of EST-contigs from Stevia were conserved to Arabidopsis thaliana followed by Oryza sativa japonica as compared to others (Fig. 3). The assembled EST-contigs (529) were further use for gene ontology (GO) based functional characterization. There were 2,361 gene ontology terms recovered. The results of GO showed 5 GO terms per contig; 76 terms, one contig (maximum) and one term, 6 contigs (minimum) were obtained in this study (Fig. 4) by analyses of 529 translated contigs though omics box. To predict the exact function of these assembled contigs sequences, contigs are grouped as per the gene ontology terminology of BP (Biological process) MP (Molecular function), and CC (Cellular component) though omics box at level 3 scheme (Fig. 5a, b and c).

Fig. 3

Assembled EST-contigs of Stevia rebaudiana in term of species distributions.

Fig. 4

EST-contigs vs. numbers of GO terms distribution obtained from Stevia rebaudiana.

Fig. 5

Gene ontology classification of EST-contig sequences of Stevia rebaudiana (a): Distribution of GO terms in the biological process category. (b): Distribution of GO terms of in the cellular process category. (c): Distribution of GO terms of in the molecular function category.

To predict the enzymes pathways in term of understanding the sophisticated functions and conveniences of the biological organization SBH method was applied for KEGG study. 884 annotated contigs were categorized by KEGG analyses resulted into six different functional groups (Table S3). 390 assembled contigs were showing signal for carbohydrate metabolism (33.07%), energy metabolism (27.43%), amino acid metabolism (14.61%), biosynthesis of other secondary metabolites (6.14%), lipid metabolism (3.07%), metabolism of terpenoids and polyketides (2.82%) and other sub-categories (Table S3). Moreover, 55 contigs belongs to cellular process containing of cell growth and death (45.45%), transport and catabolism (40.0%) and cell communication (14.55%). The 61 contigs were also involved in processing mechanism of genetic information; belongs to the 50.82% (folding, sorting and degradation), 40.98% (translation) and 6.56% (transcription) with the rest of 67 contigs play an important role in environmental processing by KEGG annotation.

The results of GO slim-plant analysis putative stress responsive genes determined by cis-regulatory element detection method is further confirmed through the STRING 10.0 indicated the presence of specific domain known for stress response (Fig. 6). The results of STRING 10.0 based on PFAM Protein Domains PF09425 (divergent CCT motif), PF06200 (tify doamin), PF,14215 (bHLH, MYC, MYB), PF00010 (helix loop DNA binding domain), PF007716(basic region leucin zipper) and PF00170 (bZIP) and IPR018467 (CO/COL/TOC1, conserved site), IPR040390 (TIFY/JAZ family), IPR010399 (Tify domain), IPR025610 (Transcription factor MYC/MYB N-terminal), IPR024097 (Basic helix-loop-helix leucine zipper transcription factor), IPR036638 (Helix-loop-helix DNA-binding domain superfamily), IPR011598 (Myc-type, basic helix-loop-helix (bHLH) domain) and IPR004827 (Basic-leucine zipper domain) were obtained (Fig. 6). The Features databases indicated 8 domains belongs to abiotic universal stress protein ath04075, ath04016 and ath00592 involved in plant hormone signal transduction, MAP kinase signaling pathway and lipid metabolism (alpha-Linolenic acid metabolism) respectively (Fig. 7a and b; Table 2, Table 3, Table 4). All identified putative abiotic stress responsive genes indicated correlation above 0.4 and mutual connection of at least 0.9. The results of co-expression analysis were confirmed using produced Protein-Protein Interaction (PPI) network by string (Fig. 8).

Fig. 6

InterProScan repeats distribution of Stevia rebaudiana contigs.

Fig. 7

Pathways predictions from Stevia contigs: A)- Plant hormone signal transduction, B)- MAP Kinase signaling pathway.

Table 2

KEG pathways obtained from Stevia EST-Contigs.

Term id	Term description	False discovery rate
ath04075	Plant hormone signal transduction	5.22E-23
ath04016	MAPK signaling pathway - plant	0.00034
ath00592	alpha-Linolenic acid metabolism	0.0074

Table 3

Protein Domain (Pfam).

Term id	Term description	False discovery rate
PF06200	tify domain	1.34E-27
PF09425	Divergent CCT motif	1.50E-23
PF00010	Helix-loop-helix DNA-binding domain	5.71E-08
PF00170	bZIP transcription factor	0.0025
PF07716	Basic region leucine zipper	0.0025
PF14215	bHLH-MYC and R2R3-MYB transcription factors N-terminal	0.0049
PF13912	C2H2-type zinc finger	0.0147

Table 4

Protein Domains and Features (InterPro).

Term id	Term description	False discovery rate
IPR010399	Tify domain	2.50E-27
IPR040390	TIFY/JAZ family	1.45E-25
IPR018467	CO/COL/TOC1, conserved site	6.05E-24
IPR036638	Helix-loop-helix DNA-binding domain superfamily	8.31E-09
IPR011598	Myc-type, basic helix-loop-helix (bHLH) domain	2.23E-07
IPR024097	Basic helix-loop-helix leucine zipper transcription factor	3.53E-05
IPR004827	Basic-leucine zipper domain	0.0032
IPR025610	Transcription factor MYC/MYB N-terminal	0.0069

Fig. 8

The interaction of the significant TFs responsive gene network developed by STRING10.0 showing gene/proteins relationship

Discussions

Stevia rebaudiana (Bertoni) one of the important medicinal plants is known for high concentrations of diterpene steviol glycosides. Stevia rebaudiana (Bertoni) leaves based ESTs (transcripts) may be a good target of diterpene synthesis. The applications of ESTs in genomics has importantly used by many investigators to identify the candidate gene (s) or transcription factor that may be increased our understanding of diterpene metabolism. In recent past, the applied genomics have been applied to understand the role of ESTs in term of stress-responsive elements predictions and its mechanism in plants system. These responsive elements of candidate gene (s) involved in stress responsive mechanisms could be assessed by comparative transcriptomics activity by using these ESTs either by RNA seq or whole genome libraries (Xu and Huang, 2018, Deokar et al., 2011). An important ESTs generated genomic approach has been implemented for stevia to identifying stress related transcription factor and genes reveled from NCBI database. The clustering of EST sequences and development of EST based contigs gave the new insight for the identifications of cis-elements or genes with their functional organizations and genome location. The 619 contigs were assembled by using the EGassembler online tools from 5,548 EST sequences. The 619 sequences were further analyzed and several classes of TFs were found (Table S1, Fig. 1). These TFs are involved in the various functions of stimulation of stress responsive genes or factor in many plants. 15 putative Transcriptional families were recognized; WD40-like (21.43%), followed at MYB-HB like, AP2-EREBP, C2H2, Hap3/NF-YB, bHLH (14.29 each) while C2C2-CO-like, CW-Zn, FHA-SMAD, Nin like, SBP3, TIFY, Tc-PD, Znf-B and bZIP (7.14% each) play an important role either direct or indirect in cellular process of plants in this study (Table S1). Several researches are already proved the role of these factor such as MYB involved in cold, salt and drought in Arabidopsis thaliana and rice, WD40-like protein involved in wheat, AP2-EREBP, in rice played an significant role under abiotic conditions (cold, drought, and salt stress) while C2H2 in rice involved in oxidative, stress, salt, cold and drought) in rice and other plants (Kong et al., 2015, Jisha et al., 2015, Xiong et al., 2014, Aprile et al., 2009, Huang et al., 2009, Varshney et al., 2009a, Varshney et al., 2009b). Similarly, other factor bZIP is involved in salt and drought stress tolerance in cotton and Arabidopsis (Kang et al., 2019, Liang et al., 2016).

In addition, 292 contigs recognized six enzyme classes hydrolases (35.49%) observed as followed by oxidoreductases, transferases, lyases isomerases ligases, and translocases with 231 sub enzymes (S2 Table). Under stress conditions, the steviol glycosides production is reduced (Fig. 2, Table S2). It is may be due to cellular pools of these enzymes the biosynthetic pathways (s) altered under adverse conditions. These results are also consonance with other findings under stress conditions (Isah, 2019, Deokar et al., 2011)

The stevia contigs through BLAST2GO program under BLASTX against the protein databases and with the help of KEGG pathway db have been implemented for prediction of genes and their functions. In-silico functional characterizations was revealed through BLASTx and the InterPro database. Total 390 contigs are play an important role in metabolism. 11 for terpenoids metabolism and polyketides, 24, other secondary metabolites biosynthesis and 129 for carbohydrate metabolism has been predicted. There are so many contigs belongs to the stress responsive genes or factors (Table S3). 67 contigs were classified into “environmental information processing.

Out of 619 some are belonging to NIR family (contig2), methyl-CpG-binding domain-containing protein (contig 46), contig 78 belongs to the heat shock protein 70 family (HSP70B), contig 215 belongs to the peroxidase family etc. (Table S2 and S4: Gene predictions and role of predicted stevia contigs). Most of these stress responsive factors and genes belongs to the abiotic stress in many plants (Zhichang et al., 2010, Du et al., 2019, Schroda et al., 1999, Hartl, 1996).

These 529 EST contigs were then subjected to GO functional classification according to the GO vocabulary of molecular function, cellular component, and biological process at Level 3 (Fig. 4). 2,361 GO terms retrieved, indicating an average of five GO terms per contig; there was a maximum of 76 GO terms for one contig and a minimum of one GO term for 6 contigs were involved in cellular functional (Fig. 5a-c).

Different correlation calculations have been used for categorizing co-expression analysis, indeed, reveals response of genes to a particular stress or more stresses with some similar consequences (Fig. 8). Furthermore, in order to confirm co-expressed genes, their co-regulation might be evaluated for having similar cis-regulatory elements by String 10.0 (Sharma, 2015, Shinozaki and Yamaguchi-Shinozaki, 2000). Therefore, to check out assumed functionally related genes with same co-regulation pattern, comparing the results with A. thaliana genome data could be useful after comparing these transcription factors obtained from contig440, contig510, contig49, contig398, contig589, contig395, contig23, contig96, contig46, contig76, contig185, contig192, contig444, contig480, contig246, contig436, contig330, contig58 and contig449 (Table 2, Table 3, Table 4) have been either belongs to direct or indirect for stress signaling. PPI network also indicated that proteins encoded by identified TFs had interactions with each other in co-expression prediction (Xu et al., 2020). In this study, the predicted stevia EST-contigs, play important role in different abiotic responsive pathways consequence of salt, drought, and temperature stress; plant signaling hormone pathways (Fig. 7, Table 2, Table 3, Table 4). The similar study have also been studied in other plants such as bHLH Transcription Factors (MYC2, MYC3, and MYC4) for jasmonate signaling in Arabidopsis (Wang et al., 2017), AP2 domain of APETALA2,DNA binding proteins play important role in the control of Arabidopsis flower and seed development (Okamuro et al., 1997), MADS box gene AGAMOUS (AG) like TFs controls reproductive organ identity in Arabidopsis (Mizukami et al., 1996) and bZIP transcription factor in Arabidopsis for hypoosmotic conditions (Tsugama et al., 2016).

Conclusions

Stresses influenced the signaling process that affect the common regulatory system and cross-talk in term of molecular changes. These molecular changes in regulatory network of plants affect the biosynthesis of metabolites and ultimately lead to stunt the plant growth. Stress tolerance in plants is generally under polygenic control and some of the genes regulating stress-tolerance responses happen to code for TFs. Therefore, TFs might be ideal candidate regulators for improving abiotic stress tolerance in crop plants. In addition to the stress-induced upregulation of the above proteins, the expression of various regulatory proteins including TFs, protein kinases, and protein phosphatases is also activated. Stevia rebaudiana is a producer of diterpenoid steviol glycoside (SGs), which is strictly regulated by plant metabolism. Due to adverse environments, plant growth is stunted and stresses influenced the signaling process that affect the common regulatory system and cross-talk in term of molecular changes resulted in less steviol glycoside synthesis in plants. These molecular changes in regulatory network of stevia plants may affect the biosynthesis of steviol glycoside metabolites and ultimately affected to yield of steviol glycoside. In this study, we have annotated stevia EST and developed the relationship of predicted contigs by genome assembly using different database and tools to predicts their role in plant signaling. After annotations we have found that the predicted TFs and their related genes are directly involved in plant metabolism and signing observed by GO analyses. We have also seen the co-expression and mutual connection of these stress responsive TFs and its related genes are involved in certain environmental factors (heat cold, light and drought). Finally, we can conclude that these factors (WD40-like protein, MYB-HB like, AP2-EREBP, C2H2, Hap3/NF-YB, bHLH (14.29% each) while C2C2-CO-like, CW-Zn, FHA-SMAD, Nin like, SBP3, TIFY, Tc-PD, Znf-B and bIP if proceeds through genetic engineering or overexpressed in stevia plants, it may lead to enhanced the plant growth and secondary metabolites accumulations.

Declaration of Competing Interest

The authors declare that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.