De novo transcriptome analysis of Justicia adhatoda reveals candidate genes involved in major biosynthetic pathway.

BACKGROUND: Justicia adhatoda is an important medicinal plant traditionally used in the Indian system of medicine and the absence of molecular-level studies in this plant hinders its wide use, hence the study was aimed to analyse the genes involved in its various pathways. METHODS AND
RESULTS: The RNA isolated was subjected to Illumina sequencing. De novo assembly was performed using TRINITY software which produced 171,064 transcripts with 55,528 genes and N50 value of 2065 bp, followed by annotation of unigenes against NCBI, KEGG and Gene ontology databases resulted in 105,572 annotated unigenes and 40,288 non-annotated unigenes. A total of 5980 unigenes were mapped to 144 biochemical pathways, including the metabolism and biosynthesis pathways. The pathway analysis revealed the major transcripts involved in the tryptophan biosynthesis with TPM values of 6.0903, 33.6854, 11.527, 1.6959, and 8.1662 for Anthranilate synthase alpha, Anthranilate synthase beta, Arogenate/Prephenate dehydratase, Chorismate synthase and Chorismate mutase, respectively. The qRT-PCR validation of the key enzymes showed up-regulation in mid mature leaf when compared to root and young leaf tissue. A total of 16,154 SSRs were identified from the leaf transcriptome of J. Adhatoda ,which could be helpful in molecular breeding.
CONCLUSIONS: The study aimed at identifying transcripts involved in the tryptophan biosynthesis pathway for its medicinal properties, as it acts as a precursor to the acridone alkaloid biosynthesis with major key enzymes and their validation. This is the first study that reports transcriptome assembly and annotation of J. adhatoda plant.

Introduction

Justicia adhatoda Linn. is a perennial shrub found in the tropical regions of Southeast Asia belonging to the Acanthaceae family. It is widely branched with pink, purple, or white flowers and can grow up to a height of 2.5 m. The plant has been widely used in the Ayurvedic and Unani systems of Indian medicine for treating various diseases such as bronchitis, common cold, asthma, cough, and tuberculosis [1]. J. adhatoda has multiple pharmacological properties that have been used to treat blood disorders, jaundice, mouth ulcers, vomiting, fever and heart-related problems [2, 3]. Phytochemical screening of J. adhatoda methanolic leaf extract revealed the presence of alkaloids, anthraquinones, cardiac glycosides, flavonoids, saponins, polyphenols, terpenoids and phytosterols [4]. The leaf extract of J. adhatoda acts as a potent anti-diabetic and anti-lipidaemic agent [5] and is also found to possess various pharmacological properties such as anti-microbial, hepatoprotective, wound healing, anti-ulcer, and anti-inflammatory activity [6]. The major alkaloids present in the leaf extracts are quinazoline alkaloids such as vasicoline, adhatodine, vasicolinone, adhavasinone, and vasicol [6]. GC-MS analysis of methanolic leaf extract of J.adhatoda revealed the presence of various bioactive compounds such as 9, 12, 15-Octadecatrienoic acid (Z, Z, Z), n-Hexadecaonoic acid, Phytol, and Amrinone phytocompounds [7]. Anisotine and vasicoline are potent inhibitors of the main protease, which is a key component for cleaving the viral polyprotein that is targeted for treating COVID-19 [8]. The development of next-generation sequencing has been useful to elucidate the secondary metabolites and its candidate genes with low labour and cost-effectiveness [9]. De novo sequencing plays a vital role in gene discovery where a reference genome is not available in a non-model organism, as it is a cost-effective method [10]. An assembly of short-read sequence data is used to identify the candidate genes involved in biosynthetic pathways and gene expression analysis [10]. Quantitative real-time PCR elucidates the level of gene expression that encodes the various components involved in the biosynthetic pathways. Despite the well-established role of J. adhatoda in the Indian medicinal system, genomic level study of the plant is not well known [11]. Thus, the present study provides a better understanding of the biosynthetic pathways and the candidate genes that are attributed to the medicinal properties of J. adhatoda.

Materials and methods

Plant material and RNA isolation

Mature and healthy leaves of J.adhatoda were collected and RNA extraction was carried out using TRIzol® Reagent (Invitrogen, USA). The total RNA extracted was treated with DNase A and purified with the RNeasy MinElute clean-up kit (Qiagen Inc., GmbH, Germany, USA). The quantity and quality of total RNA extracted were evaluated using Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, Delaware, USA). RNA integrity value was measured using Bioanalyzer 2100 (Agilent Technologies, Santa Clara, California, USA). The sample was further subjected to cDNA library preparation to carry out Illumina sequencing [11].

Library preparation and Illumina sequencing

The total RNA was made rRNA free using the Ribo-Zero rRNA removal kit (Illumina Inc., Singapore) and the remaining fraction was purified and eluted.The purified RNA sample was fragmented into short sequences using fragmentation buffer, and these fragments were used for first-strand cDNA synthesis using superscript II reverse transcriptase (Invitrogen, Carlsbad, Califonia, USA). The second strand was synthesized and purified. Poly-A tailing and adapter ligation was performed for paired-end library preparation. The primers based on adapter region were used for amplification of library for enrichment of the cDNA fragments. The quality of sequencing library was assessed using Caliper LabChip GX using HT DNA High Sensitivity Assay kit ( Caliper Life Sciences Inc., USA) Sequencing was performed with the NextSeq 500 using TruSeq v3-HS kit to generate 100 bp paired-end reads (Illumina Inc., USA).

De novo transcriptome assembly and clustering

The quality of raw paired-end reads was assessed using FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and the adapter sequences and low-quality bases were removed using Cutadapt v3.5 and Sickle v1.33 tools, respectively [12-14]. The bases with Phred score ≥ 30 reads were retained, and the de novo assembler Trinity v2.13.2 [15] was used for the transcriptome assembly, which combines three independent software modules: Inchworm assembles the data into full-length transcripts for a dominant isoform;Chrysalis clusters the contigs and constructs de Bruijn graphs for each cluster, representing the transcriptional complexity for a gene; Butterfly processes the graphs and reports full-length transcripts for alternatively spliced isoforms. CD-HIT v4.8.1 [16] was used for the clustering to remove the redundancy of assembled contigs.

Assessment of gene completeness

An online tool TRAPID was used to analyse functional and comparative de novo transcriptome data set (http://bioinformatics.psb.ugent.be/trapid_02/), which compared the unigene transcripts against PLAZA 4.5 dicots plant database with an E-value < 1E-5 for significant similarity search and annotation [17]. Considering one or more hits in the TRAPID database for full length, quasi full length, or partial length based on ORF, unigene completeness was assessed [18].

Functional annotation and gene ontology (GO) classification

The sequences of transcriptome assembly were compared against the non-redundant (nr) protein sequence database at National Centre for Biotechnology Information (NCBI) using BLASTX function from BLAST + package. The output was further analysed by Blast2GO software for retrieving the Gene Ontology (GO) terms of assembled unigenes[19]. The Kyoto encyclopedia of genes and genomes (KEGG) was performed for pathway mapping. The KEGG Automated Annotation Server (KAAS) database was also used for mapping the biosynthetic pathways of various seconday metabolites with an E-value of 1E-05. The GO terms were classified based on three different categories including Biological process, Molecular function and Cellular component.

Simple sequence repeat (SSR) detection

Simple Sequence Repeats (SSRs) of J. adhatoda were identified using MISA tool (https://webblast.ipk-gatersleben.de/misa/) with default parameters to identify mono-, di-, tri-,tetra-, penta-, and hexa-nucleotide motifs as 10, 6, 5, 5, 5, and 5 repeats respectively.

Transcript quantification

The quantification of transcripts of J. adhatoda leaf was performed using the SALMON tool [20]. It quantifies the transcripts based on the GC content, which estimates the sensitivity and abundance of differential gene expression, and quasi mapping of reads was done, where the transcripts mapped itself fast and accurately.

Validation by real-time PCR

To validate the transcriptome assembly as well as to understand the expression of unigenes from J. adhatoda, qRT-PCR analysis was performed using QuantStudio 5 Real-Time PCR system (Thermo Scientific, Wilmington, Delaware, USA) and QuantiNova SYBR Green PCR Kit (Qiagen Inc., GmbH, Germany). Control reaction without a template was included for each selected gene. The Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) from J. adhatoda was used as an internal reference gene for normalization and estimation of gene expression. A comparative Ct (2−ΔΔCt) method was used to analyse the qRT-PCR data, and fold change in gene expression was calculated using ΔCt values [21]. The experiment was repeated with three technical and two biological replicates.

Results

RNA sequencing and de novo transcriptome assembly

RNA sequencing of leaf transcriptome from J. adhatoda generated 51,639,915 raw reads and the pre-processing of raw data was performed to remove the adapter sequences, 39,134,610 high quality reads were retained with GC content of 41.69%. The raw reads were deposited at National Centre for Biotechnology Information (NCBI) Short Read Archive (SRA) database under the Accession number PRJNA84216. A total of 171,064 transcripts were assembled from raw reads with maximum and minimum unigene lengths as 15,630 bp and 201 bp, respectively with an N50 value of 2065 bp for the assembled transcriptome. The transcriptome assembly details are mentioned in Table 1 and The length distribution of unigenes is given in Supplementary Fig. 1.

Table 1

Summary of paired-end sequencing and de novo assembly of J. adhatoda transcriptome

Particulars	Numbers
Number of raw reads	51,639,915
Number of clean reads	39,134,610
No. of bases (after processing)	33,247,241
Mean Phred score	37.45
Total transcripts	171,064
Percentage of successful assembly from raw reads	64.38
Average length (bp)	1383.89
Median contig length	1096
GC %	41.69
Contig N50 (bp)	2065

Gene completeness analysis of J. adhatoda resulted in meta annotation of 30,857 (28.9%) full-length unigenes, 14,071 (13.2%) quasi full-length unigenes, 22,737 (21.3%) of partial unigenes, and 39,221 (36%) unigenes did not match to any proteins in PLAZA 4.5 dicots plant database.

Functional annotation of unigenes

The de novo assembled unigenes of J. adhatoda were annotated for sequence similarity search using BLASTX against a non-redundant protein database at NCBI with a 1E-5. The results showed 105,572 annotated unigenes and 40,288 non-annotated unigenes, of which 4772 unigenes were predicted due to inadequate genomic information in public databases. The unigene similarity search reveals that the plant has the best similarity with Erythranthe guttata (55%), Coffea canephora (6.3%), Utricularia gibba (3.4%), among others. Thus, the results indicated that J .adhatoda is more closely associated with Erythranthe guttata.

Gene Ontology (GO) classification of unigenes

GO classification was performed using the BLAST2GO tool in order to classify the genes based on gene annotation into three different categories: molecular function, cellular component, and biological process and 47 sub-categories. A total of 143,277 genes were assigned to 20 classes in the biological function, 26,732 genes in organic substance metabolic process, 26,713 genes in primary metabolic process, 21,100 genes in nitrogen compound metabolic process, 19,714 genes in the cellular metabolic process, and 11,476 genes in the biosynthetic process. In molecular function, 93,438 genes were classified into 14 sub-categories, which include the heterocyclic compound of 21,276 genes, organic cyclic compound binding with 21,276 genes, 16,105 genes with transferase activity, 13,563 genes with small-molecule binding and 12,904 genes with hydrolase activity. Additionally, 95,725 genes were assigned to cellular component category with 13 sub-categories where intracellular had the highest number of unigenes (23,055), the intracellular part 21,946 genes, 18,296 genes of intracellular organelle, and 16,066 genes of the membrane-bound organelle. (Supplementary Fig. 2).

Biological pathway analysis

The biochemical pathways of J. adhatoda were identified by mapping the unigenes to the KEGG pathway database with the help of the KAAS and BLAST2GO software with 5981 unigenes annotated to 144 biochemical pathways. A total of 6409 unigenes were assigned to the metabolic pathway, including nucleotide metabolism, carbohydrate metabolism, and amino acid metabolism having 2539, 347, and 477 unigenes, respectively (Supplementry Fig. 3). The KEGG pathway analysis of secondary metabolite biosynthesis was divided into 12 sub-categories, where the highest number of unigenes were found in sesquiterpenoid and triterpenoid biosynthesis (57 unigenes), followed by ubiquinone and another terpenoid-quinone biosynthesis (26 unigenes), terpenoid backbone biosynthesis (22 unigenes), and phenylpropanoid biosynthesis (11 unigenes) (Fig. 1).

Fig. 1

KEGG pathway unigene assignments based on secondary metabolite biosynthesis

Tryptophan biosynthesis

The major genes involved in the tryptophan biosynthesis pathway from the KEGG database are presented in Supplemnetry Table 1 and the genes identified in this pathway from transcriptome data have been depicted in Table 2. KEGG pathway map representing the biosynthesis of Tryptophan has been depicted in Fig. 2.

Table 2

Major genes involved in tryptophan biosynthesis pathway from the J.adhatoda leaf transcriptome

Gene name	E C Number	Unigene ID		Unigene length	TPM value	No. Of Unigene
Anthranilte synthase alpha	4.1.3.27	TRINITY_DN13973_c0_g1_i2	524		6.0903	1
Anthranilate synthase beta	4.1.3.27	TRINITY_DN9855_c0_g1_i4	2058		33.6854	1
Arogenate/ Prephenate dehydratase	4.1.1.48	TRINITY_DN8030_c0_g1_i1 TRINITY_DN11169_c0_g1_i1 TRINITY_DN11169_c0_g1_i2 TRINITY_DN11169_c0_g1_i3 TRINITY_DN9715_c0_g1_i2 TRINITY_DN9715_c0_g1_i3 TRINITY_DN13484_c0_g1_i1	1605 2648 3676 3139 1765 666 769		6.9556 0.8030 2.3483 8.1132 11.527 6.4298 5.9688	7
Chorismate synthase	4.2.3.5	TRINITY_DN12680_c0_g1_i3 TRINITY_DN12680_c0_g1_i5	2687 1132		1.6959 0.8975	2
Chorismate mutase	5.4.99.5	TRINITY_DN11294_c0_g1_i1 TRINITY_DN11294_c0_g1_i5 TRINITY_DN11294_c0_g1_i6 TRINITY_DN11294_c0_g1_i7 TRINITY_DN11295_c0_g1_i1	2184 1293 2287 2016 1170		8.1662 2.2369 2.1131 1.0050 4.4995	5

Fig. 2

Tryptophan biosynthesis pathway based on KEGG mapping represents the various enzymes identified from transcriptome data concerning Justicia adhatoda in each color along with its enzyme code (EC number). The KEGG pathway map is adapted from map00400 https://www.kegg.jp/keggbin/highlight_pathway?scale=1.0&map=map00400&keyword=tryptophan

Identification of simple sequence repeats (SSRs)

A total of 25,978 SSRs were identified from 106,886 sequences with 4374 sequences having more than 1 SSR, and 1951 SSRs were present in compound formation. The number of SSR loci with 8530 di-nucleotide repeats, 6802 tri-nucleotide repeats, 665 tetra-nucleotide repeats, 96 penta nucleotide repeats, and 61 hexa nucleotide repeats are represented in Supplementry Table 2. SSRs with five tandem repeats (3984) were the most common in J. adhatoda, followed by six tandem repeats (3902), seven tandem repeats (2289), nine tandem repeats (2087), eight tandem repeats (1622), and ten tandem repeats (589). Among di-nucleotide repeats, AT/AT was found to be the highest with 2968 repeats, followed by AG/CT with 2294 repeats, and tri-nucleotide repeats AAG/CTT has the highest frequency of 1737 and AAT/ATT with 1669 repeats and other motifs distributed uniformly (Supplementry Fig. 4).

The expression levels of de novo assembled unigenes of the J.adhatoda leaf transcriptome were calculated based on TPM values with the help of the SALMON tool. In the tryptophan biosynthesis pathway, the TPM values of key enzymes involved are 6.0903, 33.6854, 11.527, 1.6959, and 8.1662 for Anthranilate synthase alpha, Anthranilate synthase beta, Arogenate/ Prephenate dehydratase, Chorismate synthase, and Chorismate mutase, respectively (Table 2). The top 10 most abundant unigenes in the J.adhatoda leaf transcriptome are represented in Supplementry Table 3.

The qRT-PCR analysis was performed to analyse the expression pattern of the selected tryptophan biosynthesis genes and to validate the transcriptome assembly. The genes selected are Anthranilate synthase alpha (EC: 4.1.3.27), Anthranilate synthase beta (EC: 4.1.3.27), Arogenate/Prephenate dehydratase (EC: 4.1.1.48), Chorismate synthase (EC: 4.2.3.5) and Chorismate mutase (EC: 5.4.99.5). Anthranilate synthase alpha, Anthranilate synthase beta, Chorismate synthase and Chorismate mutase showed significant up-regulation in mature leaf when compared to young leaf and root. Arogenate/prephenate dehydratase was down-regulated in mature leaf tissue. GAPDH was used as a housekeeping gene for the gene expression analysis. Anthranilate synthase beta and Arogenate/prephenate dehydratase though with low TPM values were chosen for qRT-PCR analysis, as they were also involved in the tryptophan biosynthesis pathway (Fig. 3).

Fig. 3

qRT-PCR validation of Anthranilate synthase alpha 1 (ASA1), Anthranilate synthase beta (ASB), Arogenate/prephenate dehydratase (ANPRT), Chorismate synthase and Chorismate mutase in a young leaf (YL), mid-mature leaf (MDL), and root (ML) of Justicia adhatoda

Discussion

De novo transcriptome assembly of J. adhatoda revealed various enzymes involved in the secondary metabolite biosynthesis pathway that play a vital role in exhibiting various pharmacological properties. These pathways were identified using KEGG database where 144 biochemical pathways were identified using the assembled transcripts. Tryptophan is an essential amino acid, and a substrate, required for synthesis of serotonin, a neurotransmitter produced in the brain. The major genes identified in tryptophan biosynthesis were Anthranilate synthase alpha (EC: 4.1.3.27) and Anthranilate synthase beta (EC: 4.1.3.27). They are involved in the conversion of Chorismate into Anthranilate; Arogenate/prephenate dehydratase (EC: 4.1.1.48) converts 1-(2-Carboxyphenylamino)-1’-deoxy-D-ribulose 5-phosphate into (3-Indoyl)-glycerolphosphate, Chorismate synthase (EC:4.2.3.5) converts Chorismate into 5-O-(1-Carboxyvinyl)-3-phosphoshikimate, and Chorismate mutase (EC: 5.4.99.5) converts Prephenate into Chorismate. Anthranilate synthase is a rate-limiting enzyme involved in the terpenoid indole alkaloid pathway (TIA) [22]. It is a holoenzyme with two heterotetramers consisting of two alpha and beta subunits. The alpha subunits catalyse the formation of Chorismate to Anthranilate, which plays a vital role in the TIA pathway along with the binding site of tryptophan involved in feedback inhibition. The aminotransferases activity of the beta subunit plays a vital role in the transfer of amino groups from Glutamine to the alpha subunit of Anthrnilate synthase [23, 24]. Chorismate is the final compound of the shikimate pathway that is formed by the enzyme chorismate synthase. This enzyme catalyses the trans-1,4 elimination of phosphate from 5-enolpyruvylshikiukimate 3-phosphate, and it is of two types: fungal-type bifunctional chorismate synthase with NADPH-dependent flavin reductase and bacterial and plant-type monofunctional chorismate synthase [25, 26]. Prephenate aminotransferase and arogenate dehydratase catalyses the final step of phenylalanine production [27-29]. Anthranilate synthase is a rate-limiting enzyme that catalyses chorismate to anthranilate in the Indole pathway biosynthesis [22, 23]. Anthranilate synthase consists of two subunits in its holoenzyme such as alpha and beta subunits. Alpha subunit plays a vital role in catalysing chorismate to anthranilate by feedback inhibition. The aminotransferase activity of the beta subunit of anthranilate synthasetransfers the amino group from glutamine to the alpha subunit [24]. Overexpression of anthranilate synthase increased the level of tryptophan that acts as a precursor for various biosynthetic pathways [30]. The production of vasicinone and vasicine was enhanced by stimulating anthranilate synthase activity by increasing the tryptophan and sorbitol in the culture [31]. Chorismate mutase is a key enzyme that catalyses the formation of prephenate from chorismate for the biosynthesis of aromatic amino acids from the shikimate pathway [32]. In a study, Chorismate mutase was found to be a putative enzyme that is bifunctional and involved in the biosynthesis of phenylalanine, osmotic, and antibiotic tolerance [33, 34]. Acridone is a major alkaloidal compound present in J. adhatoda, which acts as an anti-cancer agent [35]. It is aheterocyclic alkaloid that contains tri-cyclic rings with a carbonyl group at the ninth position and nitrogen at the tenth position [36]. There are various acridone derivatives such as glyforine, thioacridones, and acronycine, all of which exhibit pharmacological properties such as antimicrobial, antipsoriatic, anti-malarial, and anti-cancer activity [37-40]. 25,978 Simple Sequence Repeats (SSRs) were identified. Many PCR based markers like RAPDs and AFLPs were usually being developed for genetic variation studies but Microsatellites with unique DNA sequences and short repeatative flanking traits can be used for polymorphism studies. The traditional techniques are quite laborious compared to the latest technology. [41, 42]. Transcripts were quantified using Salmon tool and their Transcript per million values was obtained and alternatively RSEM tool can also be used to quantify the transcripts.

Conclusion

This is the first report on de novo transcriptome assembly and annotation of the J. adhatoda and the identified genes involved in the tryptophan biosynthesis pathway related to medicinal properties. The tryptophan biosynthesis acts as a precursor for acridone alkaloid biosynthesis. These results serve as an important resource for the molecular studies on biosynthesis of various medicinally important compounds from this plant.

Below is the link to the electronic supplementary material.

Supplementary Material 1