Novel insights into structure-function mechanism and tissue-specific expression profiling of full-length dxr gene from Cymbopogon winterianus.

The 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR; EC1.1.1.267), an NADPH-dependent reductase, plays a pivotal role in the methylerythritol 4-phosphate pathway (MEP), in the conversion of 1-deoxy-d-xylulose-5-phosphate (DXP) into MEP. The sheath and leaf of citronella (Cymbopogon winterianus) accumulates large amount of terpenes and sesquiterpenes with proven medicinal value and economic uses. Thus, sequencing of full length dxr gene and its characterization seems to be a valuable resource in metabolic engineering to alter the flux of isoprenoid active ingredients in plants. In this study, full length DXR from citronella was characterized through in silico and tissue-specific expression studies to explain its structure-function mechanism, mode of cofactor recognition and differential expression. The modelled DXR has a three-domain architecture and its active site comprised of a cofactor (NADPH) binding pocket and the substrate-binding pocket. Molecular dynamics simulation studies indicated that DXR model retained most of its secondary structure during 10 ns simulation in aqueous solution. The modelled DXR superimposes well with its closest structural homolog but subtle variations in the charge distribution over the cofactor recognition site were noticed. Molecular docking study revealed critical residues aiding tight anchoring NADPH within the active pocket of DXR. Tissue-specific differential expression analysis using semi-quantitative RT-PCR and qRT-PCR in various tissues of citronella plant revealed distinct differential expression of DXR. To our knowledge, this is the first ever report on DXR from the important medicinal plant citronella and further characterization of this gene will open up better avenues for metabolic engineering of secondary metabolite pathway genes from medicinal plants in the near future.

Introduction

Citronella (Cymbopogon winterianus), an aromatic grass with large reserves of monoterpenes and sesquiterpenes, is grown for commercial and industrial purposes in tropical as well as subtropical regions of Asia, America and Africa [1]. Industrial interest in essential oils is due to their application as fragrance in perfumes, as flavor additives for use in food products or even as pharmaceutical products [2]. Citronella is traditionally known for its insect-repellent nature, especially against malaria-causing mosquito species, therapeutic properties against various diseases [3-5] and its anti-fungal property [6]. Strong antifungal activity against several species of Aspergillus, Penicillium and Eurotium due to presence of citronellal and linalool, components of citronella oil, has also been reported [7]. Citronella oil is used to expel worms or other parasites from the intestines, control muscle spasms, increase appetite, and increase urine production (as a diuretic) to relieve fluid retention [8].

The Java citronella (C. winterianus) is known for its high quality and quantity of essential oils ranging from 0.8% to 1.0% of total weight having constituents like limonene (1.8%), geraniol (∼23%), elemol (∼10%), geranyl acetate (4.0%), α-cadinol (8.0%), citronellol (∼10%), citronellal (∼35%) etc., [9]. The processors for all these isoprenoid compounds are two common 5-carbon molecules i.e., isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP). In higher plants, these two isoprenoid compounds are synthesized by two non-related biosynthetic pathways, from different precursors [10-13]. In the mevalonate (MVA) pathway, IPP and DMAPP are synthesized from mevalonic acid whereas, pyruvate and d-glyceraldehyde 3-phosphate are used as precursors for IPP and DMAPP synthesis via mevalonate-independent plastidial methylerythritol 4-phosphate (MEP) pathway [14-20].

The MEP pathway is composed of seven enzymatic steps. But the second step i.e., NADPH-dependent reduction and conversion of 1-deoxy-d-xylulose-5-phosphate (DXP) into MEP by intramolecular rearrangement is considered as most crucial [21-23] (Fig. 1). This reaction is catalyzed by 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR: EC1.1.1.267); which is the first committed step of MEP pathway [24]. It has been seen that disruption of the 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR) gene in Arabidopsis results in albino, dwarf and defects in trichome initiation and stomata closure due to insufficient supply of all the major plant hormones synthesized from IPP and DMAPP [25]. Their results reveal a critical role for the MEP biosynthetic pathway in controlling the biosynthesis of isoprenoids. DXR gene of Zymomonas mobilis has been characterized which shows highest percentage of sequence identity (i.e., 48.2%) with dxr of Escherichia coli [26]. Recently, several groups concentrated on molecular cloning and characterization of a 1-deoxy-d-xylulose 5-phosphate reductoisomerase gene from various plant species viz., Ginkgo biloba [27], Elaeis guineensis Jacq. [28], Camptotheca acuminate [29] and Salvia miltiorrhiza [30]. The expression and molecular analysis of the Arabidopsis dxr gene along with the comparative analysis of all plant DXR sequences known to date, where it possess an N-terminal transit peptide for plastids, with a conserved cleavage site, and a conserved proline-rich region at the N terminus of the mature protein, which is absent in prokaryotic DXR homologs. It has also been reported that Arabidopsis DXR participates in the control of 2-C-methyl-d-erythritol 4-phosphate pathway [31].

Fig. 1

1-Deoxy-d-xylulose-5-phosphate reductoisomerase reaction (Adopted from Takahashi et al. (1998) [23]).

The DXR is a key enzyme in MEP pathway and therefore a change in DXR activity could alter the flux of isoprenoid expressed down the pathway. Metabolic engineering of essential oil yield and composition in mint by altering expression of deoxyxylulose phosphate reductoisomerase and menthofuran synthase improves flux through the MEP pathway that leads to increased monoterpene production [32,33]. The sheath and leaf of C. winterianus also accumulates large amount of terpenes and sesquiterpenes, compound with proven medicinal value and economic uses, so sequencing of full length dxr gene and its characterization would enable prospective manipulations of the pathway that may lead to increased accumulation of active ingredients. However, the unavailability of X-ray- or nuclear magnetic resonance (NMR)-derived structures of DXR has hindered the understanding of three-dimensional architecture, mode of cofactor recognition and critical residues participate in catalysis. Therefore in the present study, we have utilized the full length dxr gene sequence (GenBank Accession number: KJ749651) from one of our current citronella leaf whole transcriptome sequencing study (BioProject ID: PRJNA263976). In this study we have used a combinatorial approach involving homology modelling, molecular dynamics simulations and docking to explain the 3-dimensional domain architecture, dynamics behaviour and mode of cofactor recognition in DXR to facilitate the understanding of catalysis of dxr in citronella. Further, tissue specific differential expression analysis of DXR was performed both by semi-quantitative RT-PCR and quantitative real time PCR in various tissues of citronella plant. The results from the present study are expected to open up better avenues for metabolic engineering of important secondary metabolite pathway genes from important medicinal plants like citronella, and also in other related medicinal plant species.

Results and discussion

Sequence analysis

The domain prediction of DXR by SMART against Pfam database revealed two putative domains viz., 1-deoxy-d-xylulose 5-phosphate reductoisomerase (Pfam: DXP_reductoisom) positioned between 81 and 209 amino acid and 1-deoxy-d-xylulose 5-phosphate reductoisomerase (Pfam: DXP_redisom_C) comprised of 223–306 amino acids. The mature 465 amino acid peptide is neutral in nature (pI of 6.44) with a molecular mass of 50.8 kDa. The aliphatic index of DXR is very high (96.5) indicating the stability of the protein over a wide range of temperature. Instability index was 35.1 (<40) DXR is classified as a stable protein. The GRAVY index of DXR is very low (0.056) indicating the possibility of its better interaction with water. To characterize the sequence of the N-terminal region of CwDXR, we aligned CwDXR with the same region of the plant DXRs known to date (Fig. 2). The plant enzyme contains an extension of 73–80 residues that is not present in the prokaryotic sequence [24]. The prediction of transit peptide for plastids with the ChloroP Program [34] in CwDXR sequence envisaged transit peptide processing site at the N terminus of a conserved Cys-Ser-X motif (49–51 position), where X means any of the hydrophobic residues Ala, Val, or Met (Met51) in CwDXR). At the N-terminal end, the sequence is weakly conserved but enriched in Ser residues, features that are typical of plastid transit peptides [35]. In contrast, the extended region at the N-terminal side (positions 52–80 of CwDXR) is highly conserved and specifically rich in Proline residues (Fig. 2). The number of Pro residues in this region ranges from 6 to 8. The consensus motif P(P/Q) PAWPG(R/T) A can be defined in the Proline rich region of plant DXR (positions 55–62 of CwDXR sequence). The primary sequence analysis suggests that CwDXR has a transit peptide for plastids, is processed at a conserved cleavage site, and contain an extended Pro-rich region at N-terminus of the mature protein. This consensus P(P/Q)PAWPG(R/T) Pro-rich motif of all plant DXR sequences may be used as a signature for DXR protein. Analysis of secondary structure elements through CONCORD revealed random coils (55.47%) dominates over helices (33.33%) followed by strands (11.20%).

Fig. 2

Multiple sequence alignment of CwDXR sequence with all the available plant DXR sequences using MEGA v6.1 and ESPript. Conserved regions are highlighted in red square boxes labelled in white whereas partially or conserved residues are marked in yellow square boxes labelled in black.

Phylogenetic analysis

To characterize the full length sequence of DXR in citronella, the polypeptide was aligned with few other homologs’ plant DXRs known to date. Multiple sequence alignment of DXR with its closest homologs revealed that initial region of the N-terminal of all the selected sequences(1–55 amino acids) are not conserved whereas, the central region and C-terminal region is strongly conserved. The 2-D phylogenetic tree inferred from NJ method revealed dichotomy with two distinct clusters (Fig. 3). Cluster-I comprised of DXR from monocots including citronella whereas, dicots formed the second cluster (Cluster-II) well supported by strong bootstrap values within their nodes. So, it can be inferred that DXR gene may have evolved due to parallel duplication from the same ancestral gene which has transferred subsequently to different lineages of plant species. The closest homology of citronella DXR has been found with Zea mays (97% identity) followed by Seteria italica (96% identity) which formed distinct cluster (Cluster I).

Fig. 3

Phylogenetic analysis of DXR from citronella with its closest homologs using Neighbor-Joining method. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the p-distance method and are in the units of the number of amino acid differences per site. Evolutionary analyses were conducted in MEGA6.

Structure–function studies

There were three top hits of DELTA-BLAST search i.e., 1R0K-A, 1Q0L-A and 3IIE-A which were considered as putative templates for DXR (Table 1). Out of the three templates, 1R0K-A (DXR of Z. mobilis [36]) showed greater sequence identity 45% with coverage 83% over other templates. Pairwise alignment of target with closest homolog (1R0K) using Multalign Program portrayed the N-terminal (1-74 residues) didn’t show any significant homology with target, so this region was not taken into consideration for further model building exercises (Fig. 4).

Table 1

Templates selected for combative modelling of DXR from citronella.

PDB ID	Source	Total score	Query coverage (%)	Sequence identity (%)	E-value	Resolution (Å)
1R0K-A	Zymomonas mobilis	388	83	45	3e-107	1.91
1Q0L-A	E. coli	406	75	45	2e-100	2.61
3IIE-A	Yersinia pestis	401	83	45	2e-100	2.21

Fig. 4

The pair-wise sequence alignment of the target DXR and template 1R0K was constructed using Multalign and ESPript. The secondary structural elements were identified from the 1R0K structure using ESPript. The α-helices, η-helices, β-sheets and strict β-turns are denoted α, η, β and TT respectively. Similar amino acids are highlighted in boxes, and completely conserved residues are indicated by white lettering on a red background.

Comparative modelling and model validation

The optimised model after a round of energy minimization was assessed for residue-by-residue stereo-chemical quality evaluation using various model validation servers. Ramachandran plot analysis using Procheck revealed that the DXR model had 93.3% Φ and Ψ angles of its residues within the core region of the Ramachandran plot as compared with 93.5% of the template, which signifies the accuracy of modelled DXR (Fig. 5). All the bond distances and angles of modelled protein lie within the allowed range of the standard dictionary values indicating that DXR model is reasonably good in geometry and stereochemistry. The Φ and Ψ distribution of Ramachandran plot of non-glycine, non-proline residues in the model and template (1R0K-A chain) are summarized in Table 2. The packing quality of each residue of the model was assessed by Verify3D Program, where the compatibility of the model residues with their environment is assessed by a score function. Residues with a score over 0.2 should be considered reliable. The score of the refined model maximally was above 0.2 which corresponds to acceptable side chain environment (Table 2). ProSA revealed a Z-score of −9.62 for modelled (Z-score of ZmDXR is −9.68) reflecting the overall quality of the DXR model. To check the degree of structural similarity of modelled target and template, we have measured the RMSD between equivalent Cα atom pairs along with a pair-wise 3-D alignment search of the template protein with the modelled structure through iPBA and MATRAS, which showed an overall identity of mere 94.70% for 375 aligned residues with a very low RMSD of 1.08 Å. The results of iPBA and MATRAS web server are in agreement with DXR and its structural homolog sharing strong structural conservation and similar in terms of structural folding. It also signifies that the generated model is reasonably good for further studies. From the above analysis, it is evident that the geometric quality of the backbone conformation, the residue contact, the residue interaction, and the energy profile of the structure are well within the limits which confirms the reliability of the modelled structure.

Fig. 5

Ramachandran plot and ProSA analysis of modelled DXR obtained using Procheck Program and ProSA Web tool.

Table 2

Comparison of Ramachandran plot statistics of DXR model with closest structural homologue 1R0K.

Ramachandran plot statistics	Modelled structure of DXR		Template (1R0K-A)
	Residues	Percentage	Residues	Percentage
Residues in most favoured regions	319	93.3	1262	93.5
Residues in additionally allowed regions	21	6.1	18	5.3
Residues in generously allowed regions	0	0.0	1	0.3
Residues in disallowed regions	2	0.6	1	0.3
Number of non-glycine and non-proline	342	100	338	100
Number of end residues (excluding Gly and Pro)	2	—	2	—
Number of glycine residues	26	—	22	—
Number of proline residues	21	—	18	—
Overall G factor	−0.03	—	0.42	—

Molecular dynamics simulation of modelled DXR

The refinement of the modelled protein was subjected to MD simulation to get an optimized and stable structure suitable for docking with cofactor. The stability and dynamic properties of DXR model was observed by MD simulation. For cross-comparison and to understand the dynamics behaviour of the modelled DXR and its closest structural homolog, we employed MD simulation for 10 ns in explicit water for both the systems. To understand the stability of both the systems, the RMSD of Cα backbone atoms, was calculated. The calculated RMSD was found to be ∼0.3 Å; whereas, the template had an average RMSD of 0.24 Å. After initial 4 ns, the system remained in the plateau state till the 10 ns. Based on the intrinsic dynamics, and improved relaxation of the model protein, the potential energy and total energy of the structure was calculated and the radius of gyration graph was plotted with respect to different time scale. As compared to the template, DXR model had a less radius of gyri. The protein remained within compact radius of gyri of 2.2 Å indicating stability nature of the protein (Fig. 6). RMSFs of Cα atoms of the model and the template showed similar type fluctuations throughout the simulation process with minute deviation at the N-terminal end. As evidenced from Fig. 6, it can be clearly observed that the (amino acids positioned in 150–300 within the catalytic domain) residues showed more deviation reflecting flexibility of this region in protein. In addition, energy profile of the both systems i.e., the total and potential energy of model as compared the template indicated the dynamic stability of the proposed model. To assess changes in the structural arrangement of elements during 10 ns simulation, each residue was minutely observed and secondary structure elements at different time scale were mapped. As evident from Fig. 6, that most of the secondary structure elements of modelled DXR remained stable which highlight the stability and reliability of the DXR model for docking studies. A representative structure was extracted from stable of portion of trajectory (last 2 ns) to study the 3-D architecture of DXR protein.

Fig. 6

The RMSD, RMSF and radius of gyration graph of the modelled CwDXR and template 1R0K: A (ZmDXR) during MD simulation. (a) RMSD of backbone Cα atoms of the CwDXR modelled structure and ZmDXR. (b) RMSF analysis of amino acid residues of CwDXR model structure and ZmDXR. (c) Radius of gyration of CwDXR modelled structure and ZmDXR during 10 ns trajectory. All the images were generated using XMGRACE software. Black and red colors represent CwDXR and ZmDXR models, respectively.

Structural study of DXR protein

Detailed structural study of DXR protein was performed to understand the arrangement of secondary structural elements in 3-D space. Modelled DXR comprised of 17 α-helices (51.15% of amino acid, 200 residues), 11 beta-strands (16.11% of amino acid, 63 residues) and other secondary structure elements (turns/coils: 32.99%, 128 residues) (Fig. 7). So as to understand the active site architecture of modelled DXR, the pair-wise structural superposition of DXR with its closest structural homolog i.e., DXR of Z. mobilis was extensively studied. The active site of modelled DXR also consists of the two pockets i.e., NADPH binding pocket and the substrate binding pocket as also reported in ZmDXR [3]. The modelled DXR bears a three domain architecture viz., N-terminal domain (residues 1–152) with 7 β-sheets and 7 α-helices, catalytic domain (residues 153–279) with 4 β-sheets and 6 α-helices and small C-terminal domain (residues 311–391) with only 4 α-helices. This domain distribution resembled the DXR from Z. mobilis, which is composed of three domains the N-terminal domain, residues 1–150, the catalytic domain comprising residues 151–299 and the C-terminal domain, residues 300–386 [32]. The N-terminal domain contains a seven-stranded parallel β-sheet, with three α-helices on each side of the sheet. The catalytic domain also belongs to the class of α/β folds; composed of a four-stranded β-sheet, and 6 α-helices and the smaller C-terminal domain consists of a bundle of four α-helices. Interestingly, the electrostatic surface potential analysis of both DXR from citronella and Z. mobilis revealed insightful variations in the charge distribution over the substrate binding site, which can be correlated with the sequence variability (variation in the key catalytic residues i.e., evident from the pair-wise structural superposition of the modelled DXR and crystal structure) and may suggest distinct substrate-binding patterns and differences in the catalytic mechanism.

Fig. 7

Homology model of DXR from citronella and its Topology. Solid ribbon representation of the DXR model coloured by its secondary structure elements. The secondary structure elements such as α helices, β strands and the N and C termini are labelled. The images were prepared with Discovery Studio and Profunc Web server.

Docking of co-factor NADPH

The residue conservation of the binding site and structural comparisons of NADPH-dependent DXR with known NADPH-dependent forms are crucial for predicting the cofactor specificity and the enzymatic mechanism. Molecular docking was performed to understand the recognition specificity and mode of cofactor binding to DXR using CDOCKER module of DS3.5. The consensus scoring system i.e., CDOCKER ENERGY, CDOCKER Interaction Energy, Ligscore1 Dreiding, LigScore2 Dreiding, PLP1, PLP2, Jain, PMF and PMF4 used for scoring various poses of the cofactor (NADPH) with the DXR is shown in Table S2. The pose with highest CDOCKER energy was selected as the best pose and it was observed that cofactor firmly binds to active site residues of the DXR protein through hydrogen bonds and hydrophobic interactions (Fig. 8). The pose with the highest CDOCKER energy was selected to represent the studied DXR ligand inside the binding cavity and employed MMFF to minimize the final structure. The NADPH is bound in an extended conformation, through a strong network of hydrogen bonds and hydrophobic contacts (Table 3), as well as a covalent bond via Asp153 residue. The residues involved in DXR–NADPH interactions are listed in Table 3. The H-bond distances were in the range of 2.37–3.35 Å. The phosphate molecule of NADPH forms two H-bonds with Asn41 (α2 helix) with an atomic distance of 3.33 Å and 2.37 Å. Further, it also formed two additional H-bonds with Nitrogen atom of Gly39 and Ser40. Again the N atom of Ser15 (α1 helix) forms two H-bonds with O2A and O2 N atoms of NADPH. Other than these, Thr13, Gly14, Ile16, Val105, Lys128, Asp153 residues forms H-bond with different atoms of NADPH molecule with an average of atomic distance of ∼2.37 Å. In addition, a minute observation revealed that when most of the residues which participated in H-bond formation with the cofactor NADPH are discrete amino acids of α-helices. Thus it can conjecture that α-helices plays a vital role in recognition of cofactor with the DXR holo-enzyme. Comparison of the modelled DXR with closest structural homolog showed that Asn41 (NADPH binding residue in template corresponds to Asn41 of ZmDXR is only conserved whereas rest amino acid residues are not conserved indicating different mode of cofactor recognition in citronella. Apart from tight H-bond network the NADPH molecule is tightly anchored in the active pocket of DXR with widespread hydrophobic contacts where most of the residues are from various β-sheets. Moreover, theoretical calculation showed that free energy of binding score (negative CDOCKER energy) of NADPH to DXR is −1043.27 kCal/mol affirms the affinity of DXR model towards NADPH.

Fig. 8

The docked complex of DXR–NADPH. The figure shows the intermolecular H-bonds and hydrophobic interactions formed between NADPH and DXR amino acids in the NADPH binding pocket that interacted with the cofactor NADPH were shown in 3-D graphical representation. The H-bonds are presented with green dotted lines with atomic distances (Å) with interacting residues. The image was prepared using Ligplus+ software.

Table 3

Residues involved in hydrogen bonding and hydrophobic interaction of NDP with DXR obtained using LigPlus+ software.

Properties	No. of hydrogen bonds	Interacting residues forming hydrogen bonds	Hydrophobic interaction forming residues
Modelled DXR	12	Thr13,Gly14, Ser15, Ile16, Gly39, Ser40, Asn41, Val105, Lys128, Asp153	Gly11, Ala38, Arg60, Gly103, Ile104, Ala108, Ala126,Asn127, Met207, Ile211, Met269

Tissue specific expression study

A distinct differential expression pattern was observed in the semi quantitative RT-PCR analysis in which DXR expression level was found to be highest in leaf sheath followed by leaf (Fig. 9), establishing the fact that DXR is a chloroplast and other plastids specific enzyme which perfectly corroborate with the earlier findings of Lichtenthaler (1999), Lichtenthaler et al. (1997), Rodriguez-Concepcion and Boront (2002) [17,15,37]. As evidence from Fig. 10, it can be clearly observed that maximum expression was seen leaf sheath followed by leaf. Although, in number of previous studies on expression profiling of dxr gene in different plants species it was reported that higher level of expression was found in leaf tissues than other tissues tested [38]. But, it can be suggested that high level of mRNA is not enough to reflect the level of enzyme activity in a particular cell type, as mRNA or protein itself is transportable to other site of action. Our result might be an indication that the expression of dxr gene is high in plastids abundant in phloem and other cells of sheath tissue. Similar results have also been obtained in a northern blot and immunolocalization experiment conducted in Catharanthus roseus for DXR expression [39]. In C. roseus the specific expression of four MEP pathway genes including DXR in internal phloem tissues was observed which was unanticipated as photosynthetic pigments are mostly localised in mesophyll tissues. But, the data obtained suggested that internal phloem parenchyma was the main site of isoprenoid biosynthesis and implies the intercellular translocation of metabolites. Recent results suggest that phytol biosynthesis in leaf might involve translocation of an isoprenoid precursor from the vascular tissue in Arabidopsis [40]. Internal phloem plays a key role in formation of precursors and decoration of monoterpenoids. Thus, several mechanisms of transport must exist to deliver isoprenoid intermediates out of the phloem parenchyma plastids and to other compartments and/or cells where they will be further metabolised or stocked [39].

Fig. 9

Semi-quantitative RT-PCR analysis DXR gene in four tissues of citronella plant.

Fig. 10

Comparative expression of DXR gene in three different tissues of mature citronella plant using ΔΔCT method of real time PCR (RQ: Relative Quantification).

Conclusion

Citronella is a highly useful medicinal aromatic plant and so increasing its level of secondary metabolite synthesis is of interest to commercial producers of its essential oils. However, the up-regulation of the key enzyme of MEP pathway i.e., DXR will be a significant step for increased essential oil production. This study is the first ever report on characterization of full length DXR in this important medicinal and aromatic plant C. winterianus. This study is purely based on in silico methods, and so combinatorial methods involving biochemical, enzyme kinetics and site directed mutagenesis will be needed to confirm these predictions. However this is a promising first step towards opening up better avenues for metabolic engineering of secondary metabolite pathways.

Materials and methods

Primary sequence analysis

The primary protein sequence (GenBank Accession number: KJ749651) was subjected to various tools viz., InterProScan, SMART (Simple Modular Architecture Research Tool) (http://smart.embl-heidelberg.de), CDD (http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) and Pfam database (http://pfam.sanger.ac.uk/) to deduce the protein family DXR gene and to explore the domain arrangement within the protein. In addition, the Protparam tool of ExPaSy proteomic server to analyse the primary structure of the amino acid sequence of DXR protein. Consensus meta-server CONCORD (http://helios.princeton.edu/CONCORD/) was used to assign secondary structure elements from its primary amino acid sequence.

Phylogenetic study

Primary amino acid sequence of DXR was subjected to BLAST search against non-redundant (NR) database of NCBI. Sequence producing significant alignment (cut-off identity of ⩾85% and E-value of 0) with the DXR of citronella were aligned using Multalign and displayed using ESPript. Finally, 2-dimensional phylogenetic tree was established implementing Dayhoff model of substitution using Neighbor-Joining method in MEGA v6.1.0 with a bootstrap value of 1000 iterations.

Theoretical modelling of DXR

Homology model of DXR protein was generated using MODELLER v9.12 based on suitable template 1R0K-A which was chosen by using fold-recognition server’s viz., 3D Jury, Pcons.net, GeneSilico and Geno3D. Further to ensure the correctness in identifying templates for homology modelling through fold recognition servers, DELTA-BLAST (Domain Enhanced Lookup Time Accelerated BLAST) was used against protein data bank (PDB). Based on the optimised target-template alignment, modeller facilitated in the development 200 models, of which, the model with the lowest discrete optimized protein energy (DOPE) score and low RMSD (after aligning them on corresponding Cα atom pairs of the template) was selected for further refinements (i.e., loop and side chain refinement) using Discovery Studio3.5 (Accelrys, Inc., San Diego, USA). The refined model was validated by Procheck tool embedded in SAVES server to quantifies the amino acid residues in the available zones of Ramachandran plot to assess the stereo-chemical quality of the model. Energy profile of modelled DXR was analysed using ProSA tool. The root mean square deviation (RMSD) of the corresponding Cα-atom pairs of both the model and template was calculated using iPBA web server.

Molecular dynamics simulations

For further refinement the model, modelled DXR was subjected to molecular dynamics (MD) simulations using GROMOS96 43A1 force field in GROMACS 4.6 package [41-43]. For cross validation, we have used two systems i.e., modelled DXR as well its closest structural homolog ZmDXR (template: PDB ID: 1R0K). The systems were solvated with SPC3 water molecules and neutralized by adding equal number of counter ions (Na+/Cl−) with a concentration of 0.15 mol. Both the models were placed in an octahedron box maintaining a distance of 1.2 nm from all the directions of protein and periodic boundary. Then the solvated systems were subjected to energy minimization using steepest descent algorithm (2000 steps max.) followed by a 100 ps position restrained MD simulation at constant pressure (1 atm) and temperature (300 K). The long-range electrostatic interactions were computed using the particle-mesh Ewald (PME) [44] with a cut-off value of 0.9 nm for van der Waals and 1.4 nm for electrostatic interactions. All bond lengths were constrained using LINCS algorithm. The trajectories were analysed using GROMACS analysis tools and the structures were analysed using DS3.5 and PyMOL (www.pymol.org).

Molecular docking of NADPH using CDOCKER

The average coordinate of DXR protein was extracted from stable portion of the trajectory and cofactor NADPH were prepared using CDOCKER module following prepare protocol in DS3.5 prior docking. For docking we adopted the methodology from an earlier study by Dehury et al. (2014) [45]. The ligand was docked within the active site of enzymes and receptors by implementing CHARMm-based molecular dynamics method [46,47] which (binding pocket) was obtained by pair-wise alignment of DXR with the template. Based on the consensus scoring system of CDOCKER module, best ten poses were selected and binding energy of the best 10 poses of NADPH was obtained using “Binding energy module” of DS3.5. The pose with lowest energy, and more number of interacting H-bonds was selected for displaying the enzyme-cofactor interaction using LigPlot+ software.

Tissue specific expression studies

The differential expression of the DXR gene in three different tissues of citronella (leaf sheath, leaf and root) was tested by semi-quantitative RT PCR and quantitative real time PCR. The primers were designed by primer3 (http://primer3.ut.ee/) with GC% of 60–65% and melting temperature ranging from 60–65 °C. The first strand cDNA was synthesized by using PrimeScript™ 1st strand cDNA synthesis kit (Takara, Clontech) as per manufacturer’s instruction. The 10 μl reaction mixtures with 125 ng first strand cDNA with 1X Taq polymerase buffer (Genei, India), 1U Taq polymerase buffer (Genei, India) 10 pmol both forward and reverse primers (Sigma) and 5 mmol dNTP mixure (Invitrogen) were subjected to semi-quantitative RT PCR in GeneAmp Thermo Cycler (Applied Biosystem, USA). The amplification condition were 5 min at 94 °C, followed by 35 cycles of 1 min at 94 °C 1 min at 67.5 °C 1 min at 72 °C and final extension was done for 10 min at 72 °C. The amplified products were visualized in 2% agarose gel stained with EtBr. In this study rice GAPDH gene was used as reference. The relative expression of DXR gene in the same 3 tissues were checked with quantitative real time PCR by ΔΔCT method [48] on StepOnePlus Real-Time PCR System (Applied Biosystem) using the SuperScript III Platinum SYBR Green One Step qRT-PCR with ROX Kit (Invitrogen, CA, USA) according to manufacturer’s instruction. PCR amplification was performed under the following conditions: 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min finally melting at 95 °C for 15, 60 °C for 1 min. The gene expressions were normalized against an internal reference gene, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and root tissue was arbitrarily chosen to be the calibrator of tissue gene expression.