Multi-tissue transcriptome analysis using hybrid-seq revealed potential genes and biological pathways associated with azadirachtin A biosynthesis in neem (Azadirachtin indica)

doi:10.21203/rs.2.23446/v2

Download PDF

Research article

Multi-tissue transcriptome analysis using hybrid-seq revealed potential genes and biological pathways associated with azadirachtin A biosynthesis in neem (Azadirachtin indica)

https://doi.org/10.21203/rs.2.23446/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 28 Oct, 2020

Read the published version in BMC Genomics →

You are reading this older preprint version

Read the latest preprint version →

Background: Azadirachtin A is a triterpenoid from neem tree exhibiting excellent activities against over 600 insect species in agriculture. The manufacture of azadirachtin A depends on extraction from neem tissues, which is not ecofriendly and sustainable. The low yield and discontinuous supply impeded the further application. The biosynthetic pathway of azadirachtin A is still unknown.

Results: We attempted to explore azadirachtin A biosynthetic pathway and identified key involved genes by analyzing transcriptome data of five neem tissues through hybrid-seq (Illumina HiSeq and Pacific Biosciences Single Molecule Real Time (PacBio SMRT)) approach. Candidates were firstly screened by comparing expression level within five tissues. After phylogenetic analysis, domain prediction and molecular docking, 22 candidates encoding 2,3-oxidosqualene cyclase (OSC), alcohol dehydrogenase (ADH), cytochrome P450 (CYP450), acyltransferase (ACT) and esterase (EST) were proposed to be potential genes involved in azadirachtin A biosynthesis. Among them, two unigenes encoding homolog of MaOCS1 and MaCYP71CD2 were identified. An unigene encoding complete homolog of MaCYP71BQ5 was first reported. Accuracy of the assemblies were verified by Quantitative Real-Time PCR (qRT-PCR) and full-length cloning PCR.

Conclusions: By integrating and analysis transcriptome data from hybrid-seq technology, 22 DEGs were finally selected as candidates involved in azadirachtin A pathway. The obtained reliable and accurate sequencing data provided important novel information for understanding neem genome. Our data shed new light on the understanding of other triterpenoids biosynthesis in neem trees and provide reference for exploring other valuable natural product biosynthesis in plants.

Epigenetics & Genomics

Plant Molecular Biology and Genetics

azadirachtin A

natural insecticides

secondary metabolism

triterpenoid biosynthesis

transcriptome

neem

With the increasing concern on chemical pesticides threat to crop protection program globally, more attention is being paid towards bioactive, biodegradable plant or microbial based biopesticide. Azadirachtin A, the major insecticidal ingredient in neem-based products, exhibits excellent activities against over 600 kinds of insect species [1] in agricultural area [2]. It is processed by insects as natural hormone and introduces antifeedant, repellent and growth inhibiting behavior [3]. Azadirachtin A-based pesticides are biodegradable, environmentally friendly and non-toxic to mammals, plants and birds. Due to these general superiorities, the agriculture segment accounted for the highest share (40%) in total revenues of neem extract products market. It is expected to grow from 653 million dollars in 2015 to 1.8 billion dollars in 2022 with high annual growth rate of 16.3% [4]. However, current supply of azadirachtin A mostly depends on extraction of neem seed [5] and leaf. Limited by the distribution of neem materials, the production of azadirachtin A through complex and low-yield extracting approach is far from meeting the marketing demand.

As a potential insecticide, its synthesis had been investigated in the last few decades. However, its molecular architecture [6] complexity was the main obstacle that held back the advances in azadirachtin A biosynthesis. After 20 years investigation, the complete chemical synthesis of azadirachtin A was finally accomplished [7] in 2007. While the low productivity (0.00015%) of the 71 steps-reaction still limited the applications in industry.

Being a triterpenoid from neem, azadirachtin A is derived from 2,3-oxidosquanele. The downstream biosynthesis pathway (pathway from the biosynthesis of scaffold to azadirachtin A) is still unclear after decades investigation. One attempt about its biosynthesis investigation was the feeding experiments in 1971. [³H] euphol and [³H] tirucallol were used to incubate with neem leaf crude. Euphol rather than tirucallol was more efficiently incorporated into nimbolide (a limonoid with structural similarity to azadirachtin A) [8]. However, no enzyme producing tirucallol or euphol was isolated or characterized from neem.

Up to now, several resources of gene information of neem, such as cDNA library [9], expressed sequence tag library[10], draft genome[11, 12] and transcriptome data[13, 14, 15] have been available in the public database. The public has access to these resources in the database of NCBI and EBI according to the accession numbers. Genes involved in MVA, MEP pathway [14, 16] or 2,3-oxidosquanlene biosynthesis [11, 15, 16, 17] in neem were identified. Several specific genes in neem were found after comparative genome analysis of neem with Arabidopsis thaliana, Oryza sativa and Citrus sinensis [11]. Through comparing gene expression within different neem tissues, genes encoding CYP450s were identified from fruit [16] and leaf [17]. However, none of the isolated genes had been functionally characterized. It was not until the key enzyme that converting 2,3-oxidosqualene into tirucalla-7,24-dien-3β-ol, AiOSC1 was isolated from neem transcriptome [18]. This was the first functional characterized AiOSC1 reported to be related to azadirachtin scaffold biosynthesis. Two full-length CYP450s subsequently catalyzing tirucalla-7,24-dien-3β-ol into melianol were also isolated from C. sinensis var. Valencia and Melia azedarach. Homologs of MaCYP71CD2 in neem (AiCYP71CD2) was identified, while homolog of MaCYP71BQ5 (one characterized CYP450 from M. azadarach) [18] in neem was fragmented.

Although the genome and transcriptome data were easy to download from SRA database, the submitted files were all raw data without assembly or annotation. This caused difficulties for people who are not familiar with the processing of high-throughput genome or transcriptome data. Besides, samples for transcriptome or genome analysis in previous report were growing in India. Genetic information within same species in different living area or developmental stage were different and this caused the differences in gene expression. Hens, five neem tissues (fruit, leaf, stem, flower and root) were sampled for transcriptome sequencing. Among five tissues, fruit with the green hard seed was reported to contain the highest azadirachtin A in all the developing stage of fruit [14, 19]. Leaf from neem tree in China contained azadirachtin A at the value of 969.9 μg g⁻¹[20]. The rank of azadirachtin content in different tissues was consistent with that in previous report (seed kernels, 0.03%; leaves, 0.9x10^-3%; bark, 0.5x10^-3%; root, 0.3x10^-3%; stem 0.2 x10^-3%) [21]. They were chosen as the higher azadirachtin A group and used for mining synthetic genes.

With the advances in high-throughput sequencing technology, the third-generation sequencing represented by Biosciences Single Molecule Real Time (PacBio SMRT [22]) and Oxford Nanopore was applied in academic researches. Owing to its high error rate of sequencing in its longer reads (15%) as well as the accuracy of reads from Illumina, a new method called hybrid-seq [23] was generated and it bring together the best of two sequencing technologies. The obtained longer reads was corrected by short but accurate reads from Illumina. As reported by Koren [24], only 41 (0.1%) of the PacBio reads exactly matched the annotated exon structure before correction in genome assembly. This number will rise to 12, 065 (24.1%) after correction with short reads. Error correction method was used to produce de novo assembly of the Saccharomyces cerevisiae. The generated contig N50 length is more than ten times greater than an Illumina-only assembly (678 kb versus 59.9 kb), and has >99.88% consensus identity when compared to the reference genome [25]. The successful applications of hybrid-seq in genome refining, isoform identification laid strong foundations for our study.

While three steps involved in azadirachtin A biosynthesis had been figured out, the rest part of the downstream pathway was still unexplored. Based on the metabolites in Neem Metabolite Structure Database (http://www.vmsrfdatabas e.org/index.php) as well as the distribution of metabolites (Table S1) in neem tissues, a putative biosynthetic pathway for azadirachtin A (Figure 1) were proposed. Five kinds of putative enzymes (OSC, CYP450 [26], ADH [27], ACT and EST [28]) participating in were proposed based on the putative biosynthetic pathway. Candidate mining were performed through the work flow of gene mining (Figure S1). Extensive bioinformatics analysis of unigenes further provided twenty-two candidates (Table 1) including one OSC unigene, two ADH unigenes, twelve CYP450 unigenes, two ACT unigenes and five EST unigenes after phylogenetic analysis, domain prediction and molecular docking. Among them, three transcripts encoding complete AiOSC1, homologs of MaCYP71CD2 and MaCYP71BQ5[18] were also found in our study. Unigene contains complete ORF encoding homolog of MaCYP71BQ5 was first reported. Qualitative real-time PCR and full-length cloning PCR were used for verifying unigenes expression level (FPKM) and transcript sequence accuracy. The obtained candidates could be used as an important resource to investigate the catalysts responsible for essential biochemical reactions in azadirachtin A biosynthesis as well as triterpenoids metabolism in close species of neem.

Table 1. Summary of candidates involved in azadirachtin A biosynthesis in neem

Classification	Unigene	Gene accession number	Homologs
OSC	transcript/14449	gi\|443299067	AiOSC1*
ADH	transcript/18833	gi\|572153023	cinnamyl-alcohol dehydrogenase
ADH	transcript/19291	gi\|572153023	cinnamyl-alcohol dehydrogenase
CYP450	transcript/17636	gi\|225458053	CYP83B1
	transcript/17854	gi\|641826901	CYP77A3
	transcript/16057	gi\|567902124	flavonoid 3'-monooxygenase
	transcript/17284	gi\|567889747	CYP89A2
	transcript/16777	gi\|590722535	brassinosteroid-6-oxidase
	transcript/17001	gi\|590722535	brassinosteroid-6-oxidase
	transcript/16577	gi\|568834016	flavonoid 3'-monooxygenase
	transcript/16950	gi\|645239614	cytokinin trans-hydroxylase
	transcript/17057	gi\|567868115	flavonoid 3'-monooxygenase
	transcript/16971	gi\|568825869	MaCYP71CD2*
	transcript/16742	gi\|568830413	MaCYP71BQ5*
	transcript/16198	gi\|568868580	CYP51G1
ACT	transcript/17792	gi\|567873443	BAHD acyltransferase
ACT	transcript/18214	gi\|641830965	BAHD acyltransferase
EST	transcript/19188	gi\|568822600	SGNH_plant_lipase
	transcript/19882	gi\|568867507	Acetyl esterase/lipase
	transcript/19697	gi\|225440163	Acetyl esterase/lipase
	transcript/19748	gi\|641846434	Acetyl esterase/lipase
	transcript/18100	gi\|568835134	pectin acetylesterase

Homologs with * represent the functionally identified genes in the previous report.

Plant materials

All fresh and healthy tissues (root, leaf, stem, flower and fruit containing seed）from neem (Azadirachtin indica, A. indica) tree were randomly sampled from the park of Hainan University, followed by transcriptome analysis on September 26, 2018. All samples for transcriptome sequencing were harvested from three plants for RNA extraction. Tissues were gently rinsed and subsequently cut into small pieces. All materials were immediately frozen in liquid nitrogen and stored at -80℃ before use. The A. indica was authenticated by Prof. Zhiwei Wang and Dr. Yaqing Wei.

RNA extraction

Total RNA was extracted using the RNA plant Plus Reagent (TianGen, Beijing, China) according to the manufacturer's protocol. The extracted RNA concentration and integrity were assessed using the RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent, CA, USA). For PacBio sequencing, the extracted RNA concentration and integrity were assessed using the Fragment Analyzer system (Agilent, CA, USA). The A260/A280 ratio ranging from 1.9 to 2.0, concentration above 285 ng/μL, and RNA integrity number (RIN) greater than 8.0 were used for subsequent cDNA library construction.

cDNA library construction and hybrid-sequencing

A total amount of 5 μg RNA was used for cDNA library construction. For Illumina HiSeq sequencing, oligo (dT) beads were used to isolate poly(A)⁺mRNA. The paired-end libraries with an insert size of approximately 250 bp were constructed. All libraries were sequenced commercially to generate paired-end reads of average length 150 bp on the Illumina HiSeq 2000 sequencing platform (Hiseq 2000 V3) according to the manufacturer’ protocol by Beijing Genomics Institute (BGI-Shenzhen, China).

For PacBio SMRT sequencing, 1000 ng of mRNA from each tissue was pooled for cDNA library construction. Double-strand cDNA was synthesized according to SMARTer PCR cDNA Synthesis Kit (Clontech). DNAs were selected by BluePippin^TM and ranged into four sizes: 1-2, 2-3, 3-6 and 5-10 kb. DNAs after the second large-scale PCR were used as template prepared for SMRTbell library for sequencing. The amount of generated base was about 12 GB to cover all transcripts in sample.

HiSeq reads were filtered by discarding the reads with adaptor sequences, reads with ambiguous “N” bases larger than 5% and low-quality reads. and the filtered reads were then assembled using Trinity (v2.0.6) with default parameters to generate contigs. These contigs were then processed by sequencing clustering software TGICL (v2.0.6) for redundant Trinity assembled contig removal. Raw PacBio SMRT reads were processed using SMRT Analysis Server (v2.3) for full-length transcripts generation. The obtained transcripts were corrected with filtered HiSeq reads using LSC [29] software and subsequently filtered with CD-HIT-EST for the removal of redundant Trinity generated fragments. Finally, the calibrated transcripts were assembled into unique putative transcripts (including contigs and singletons) and unigenes were characterized for subsequent analysis.

Annotation and differential gene expression analysis

The unigenes were annotated based on sequence similarity using BLASTX against five databases, including non-redundant protein database (Nr), SwissProt, Clusters of Orthologous Groups of proteins (COG), and Kyoto Encyclopedia of Genes and Genomes protein database (KEGG). Pfam annotation for the unigenes was finished using the HMMER 3.0 package. Sequence description for each unigene were transferred from homologous BLAST hits with E-value<10^-5. Gene Ontology (GO) terms were assigned based on the top BLAST hit using Blast2GO. Genes were obtained by BLASTn using non-redundant nucleotide sequence database (Nt). Functional enrichment of the assigned GO terms was calculated and analyzed via WEGO software. The distribution of gene functions was illustrated by biological process, cellular component and molecular function.

Clean reads were mapped to unigenes using Bowtie2 (v2.2.5), and then gene expression level was calculated with RSEM (v1.1.12). To compare the difference of gene expression among different samples, the FPKM (Fragments per kilobase per transcript per million mapped reads) method was used [30]. DEseq2 was used to identify DEGs (absolute value of log₂ fold-change≥1) after correction of p-values (adjusted＜0.05) using Benjamini-Hochberg (false discovery rate, FDR≤0.001). Higher expressional unigenes characterized from leaf and fruit (azadirachtin A stimulating tissues) libraries were used for candidate mining.

Analysis of phylogeny, domain architecture of unigenes

The SwissProt database was queried to retrieve all reviewed sequences of alcohol dehydrogenase, CYP450, acyltransferase, and esterase. These sequences were downloaded in FASTA format and aligned with all four kinds of candidate using Clustal W algorithm with default parameters. Phylogenetic analysis based on multiple alignments of protein sequences was done using the Neighbor Joining [31] method as implemented in MEGA7 and the phylogenetic trees were visualized on iTOL (https://itol.embl.de/). Accessions of these protein sequences used in phylogenetic analysis were provided in Additional file 1. The protein sequences of the candidates were also searched against Pfam database in order to get domain architecture information complementary to that provided by SwissProt.

Molecular modelling and docking for enzyme-substrate analysis

Modelling of five CYP450 proteins encoded by unigenes was performed by using Phyre2 web portal using fold recognition method [32]. Docking analysis was performed to predict the interactions of four triterpenoids (tirucalla-7,24-dien-3β-ol, azadirone nimbin and nimbolide) as substrate for the CYP450 proteins using software AutoDock 4.0 (http://autodock.scripps.edu). For covering acting domains present in CYP450 protein, grid spacing was maintained at 0.375 Å. Genetic algorithm (GA) was applied as searching parameter with 10 GA runs and population size was set to 150, energy evaluations was set to maximum 25,00,000 with considering the maximum number of generations to 27,000. The most favorable docking conditions were in the form of lowest binding energy conformations with H-bonds in cluster. Phymol 2.3 software was used for better analysis of interaction in protein-ligand complexes obtained from AUTODOCK4 [33] software.

Validation of hybrid-seq by quantitative real-time PCR and full-length cloning PCR

Ten transcripts were selected to validate their expression from hybrid-seq by quantitative real-time PCR (qRT-PCR). Total RNA was processed with RNase-free-DNase I (TianGen, Beijing, China) following the manufacturer's instructions to eliminate potential DNA contamination. First strand cDNA was synthesized using GoScript^TM Reverse Transcription System (Promega, Canada). The reactions were performed in triplicate using 2 μL diluted cDNA template in a 20 μL total volume. qRT-PCR was performed in 96-well plates on a Bio-Rad CFX96 real-time PCR system (Bio-Rad, CA, USA) using SYBR Green Mix (Bio-Rad, CA, USA). A two-step cycling program was performed, comprising an initial 95℃ polymerase activation for 3 min, followed by 40 cycles of 95℃ for 10s, then 60℃ for 30s. The melting curve was obtained by heating the amplicon from 65℃ to 95℃ at increments of 0.5℃ per 5s. The actin gene was used as an internal control to normalize all data. The relative quantitation (ΔΔCt) method was used to evaluate differences between the tissues for each gene examined. Data analysis was performed using GraphPad Prism version 5 for Windows (GraphPad Software, Inc. La Jolla, CA, USA). The primers for qRT-PCR reactions were listed in Additional file 2.

Ten unigenes were selected to verify the accuracy of sequence obtained from hybrid-seq by full-length cloning PCR. Obtained cDNA (2 μL) was used in a 50 μL PCR reaction containing 2 μL of forward primer (10 μM), and 2 μL of reverse primer (10 μM) and 25 μL of I-5™ 2×High-Fidelity Master Mix (TsingKe Biotech, Beijing, China). The PCR product was purified using GeneJET PCR Purification Kit (Thermo, USA) and assembled into pJET vector using CloneJET PCR Clone Kit (Thermo, USA). The constructs were then transformed into Trelief 5α Chemically competent cell (TsingKe Biotech, Beijing, China) and the amplified constructs were sequenced by Genewiz Company (GENEWIZ, Beijing, China). The primers for full-length cloning PCR were listed in Additional file 2.

Hybrid-seq transcriptome analysis and data validation

To generate a comprehensive overview of neem transcriptome, total RNAs were extracted from leaves, fruits (containing seeds), roots, stems and flowers. The approach we used to obtain the transcriptome data was hybrid-seq technology which combine the Illumina Hiseq and PacBio SMRT for sequencing and correct the errors in long reads with short reads [24]. Different neem tissues were sequenced separately using Illumina Hiseq platform and generated 41.14, 41.35, 40.60, 40.59 and 40.64 million clean reads, respectively. PacBio SMRT platform produced 6.75 million clean reads. After calibration with short reads from Illumina platform, the assembled unigenes were corrected with an N50 of 5076 bp and mean length of 3607 bp (Table 2). The obtained assembly were 2.5 times and 3 times longer than in the previous report [14]. The length of these unigenes ranged from 500 to 6001bp, and the majority (over 55.5%) were distributed in 4501 bp and above (Figure 2a).

To confirm the accuracy of the hybrid-seq (FPKM) results, we selected 10 unigenes and used qRT-PCR to determine their relative expression (Figure S2). Good consistency between the qRT-PCR and FPKM were observed except for transcript/19882 and transcript/16577. Their expression from RNA-seq (FPKM) and relative expression level by qRT-PCR in root and stem were inconsistent. The inconsistency between qRT-PCR and FPKM in transcriptome analysis have also been reported by other researchers. Eight out of fifty-eight genes expression level were inconsistent with their FPKM when Zhang validated by qRT-PCR from his transcriptome [34]. In another case, when comparing gene expression fold changes between MAQCA and MAQCB samples, about 85% of the genes showed consistent results between RNA-sequencing and qRT-PCR data [35]. The inconsistency between FPKM and qRT-PCR may result from reasons below: Although realtime-PCR and RNA-seq are both used to measure gene expression, the unit of measurement [36] as well as the computing method are different for FPKM and qRT-PCR. Many factors affect the accuracy of FPKM and qRT-PCR. Bias in PCR amplification [37] and RNA-seq library preparation [38]and sequencing add noises in RNA-seq. The quality of the mRNA, amplification efficiency and the choice of reliable internal controls referred to as reference genes affect the accuracy of qRT-PCR [39]. Therefore, the consistency of transcript/16577 and transcript/19882 between qRT-PCR and FKPM was acceptable. More further examinations need to do on these two unigenes.

Ten unigenes were cloned by PCR used primers listed in Additional file 2. According to the sequencing results of the cloned unigenes, every of them was 100% identical to the sequences obtained from hybrid-seq platform. Among them, four unigenes attracted our attentions since they contained complete ORF of the genes. One of them(transcript/14554) was annotated as NADPH-cytochrome P450 reductase 2 of neem. The cloned transcript/14449 was found to encode an OSC consist of 760 amino acids. The length of two CYP450 unigenes (transcript/16971 and transcript/16742) was 1536bp and 1527bp, and they encoded protein consist of 511 and 508 amino acids, respectively. Detailed sequencing results indicated that the generated unigenes in our study were of good accuracy and therefore reliable for further analysis.

Functional annotation and classification of unigenes

In total, 19907 unigenes (98.54% of 20201 unigenes) were annotated in at least one database (Table S2). The annotated unigenes were compared to known nucleotide sequences of other plant species, which were best matched to the known nucleotide sequences from C. sinensis (52.43%), followed by Citrus clementine (23.49%), Theobroma cacao (2.44%), Vitis vinifera (2.4%) and other (19.23%).

COG and GO classification were used to further evaluate the completeness and effectiveness of the neem annotation. All 17634 unigenes (87.3% of 20201 unigenes) were classified into 25 functional Clusters of Orthologous Groups (COG) clusters (Figure 2b). Of those, 2610 unigenes (14.8% of the total 17634 classified unigenes) were categorized into general function prediction only cluster, which formed the largest group, whereas the clusters for replication, transcription, recombination and repair followed closely. Although only 378 unigenes were categorized into the “Secondary metabolites biosynthesis transport and catabolism” cluster, they may play important roles in providing precursors in secondary metabolite biosynthesis.

In total, the 74905 assembled unigenes (66.3% of 113008 unigenes obtained in Hiseq sequencing analysis) were assigned at least one of the 55 GO terms (Figure 2c). Of those, these unigenes were predominantly assigned to the metabolic process (GO:0008152) and cellular process (GO:0009987). The unigenes categorized in the molecular function category were predominantly associated with catalytic activity (GO:0003824) and binding functions (GO:0005488). The unigenes categorized in the cellular component category were predominantly associated with membrane (GO: 0016020), cell part (GO: 0044464) and cell (GO: 0005623). These findings showed that the main COG and GO classifications for the fundamental biological processes were identified.

The KEGG pathway database was used to systematically evaluate the gene biological functions participated in different pathways. A total of 16778 unigenes were matched to the database and assigned to 135 KEGG pathways (Table S3). Metabolic pathways (3590 unigenes, 21.4%) and biosynthesis of second metabolites (1560 unigenes, 9.3%) were the two dominant categories. In the biosynthesis of secondary metabolites category (Table S4) in neem, subcategories of flavonoid biosynthesis, terpenoid backbone biosynthesis (Table S5), steroid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis (Table S6), diterpenoid metabolism and carotenoid biosynthesis were included. There were 242 unigenes involved in the metabolism pathways of terpenoids and polyketides.

Differential expression analysis of unigenes in neem

In order to find candidate genes in azadirachtin A biosynthesis pathway, all transcripts had been identified and annotated and mapped in different pathways. Genes expressed higher in azadirachtin A stimulating-organ fruit and leaf would be more likely to be the involved-genes in azadirachtin A biosynthesis. The expression level of unigenes was calculated by FPKM. The number of up-regulated DEGs in leaf and fruit comparing to other tissues were 219 and 397, respectively. All these DEGs were used for mining candidate involved in azadirachtin A biosynthesis.

First screening of candidate genes involved in azadirachtin A biosynthesis

According to the putative azadirachtin A pathway in Figure 1, tirucalla-7,24-dien-3β-ol was assumed as the scaffold formed from 2,3-oxidisqualene. A few steps such as hydroxylation and furan ring formation occurred after scaffold formation. The hydroxyl groups were then either oxidized to acid or acylated or esterified to esters forming limonoid compounds like azadirone or nimbin. Azadirachtin A was finally obtained after modifications on azadirone or nimbin. ADH, CYP450, ACT and EST were supposed to be involved in azadirachtin downstream pathway and their encoding-unigenes were chosen for candidates.

As a triterpenoid, the first step of azadirachtin A biopathway was the formation of its scaffold catalyzed by OSC. Among all unigenes involved in terpenoids biosynthesis, eight detected unigenes were annotated as OSC, including transcript/1784, transcript/1866, transcript/8176, transcript/8892, transcript/9751, transcript/14584 and transcript/19700 and transcript/14449. Among them, only transcript/14449 expressed higher in fruit. While after phylogenetic analysis (Figure 3) of transcript/14449 with several characterized OSCs from other plants, transcript/14449 was grouped with AiOSC1[18], a newly characterized OSC from neem catalyzing the formation of tirucalla-7,24-dien-3β-ol, while other unigenes were grouped with cycloartenol synthase. After DNA sequence analysis with AiOSC1(Figure S3), transcript/14449 is 100% identical to AiOSC1, which means that these two genes were the same gene. It indicated that transcript/14449 could be a candidate gene for producing azadirachtin A scaffold.

Among all the DEGs in fruit and leaf transcriptome data, DEGs encoding ADH, CYP450, ACT and EST were discovered. There were sixteen DEGs encoding ADH. Four up-regulated DEGs encoding ADH were selected for further screening. As for CYP450, sixteen DEGs encoded CYP450 and ten DEGs were selected. Similarly, thirteen and nineteen DEGs encoded ACT and EST respectively and there were four and twelve DEGs up-regulating in leaf and fruit tissues, respectively. DEGs with same sequence or sequences within 150 amino acids were excluded in further screening.

Further screening of the other four enzymes through phylogenetic analysis and domain prediction

Phylogenetic trees (Figure 4) and protein domains prediction were used for further screening of these DEGs. Transcript/22186, transcript/18833 and transcript/19291 were respectively grouped with CADH4[40] and CADH1, which catalyzed biosynthesis of cinnamaldehyde from cinnamyl alcohol. Transcript/18482 was grouped with ADHX [41] which had activity to primary and secondary alcohols. Transcript/18833 and transcript/19291 contained PLN02514 domain which was also found in cinnamyl-alcohol dehydrogenase [42]. Transcript/22186 contained nsLTP2[43] domain which was contained in non-specific lipid-transfer protein. Transcript/18482 contained GxGxxG motif [44] found in S-(hydroxymethyl) glutathione dehydrogenase. Transcript/18833 and transcript/19291 were selected as candidates for further examination.

According to the phylogenetic analysis of unigenes encoding CYP450, they were divided into two clades. Among them, five transcripts (transcript/16057 and transcript/16577, transcript/16950, transcript/16777 and transcript/17001) were grouped with members in CYP71[45] and CYP72[46] clades, which members were reported to be involved in terpenoids biosynthesis. Transcript/16971 was grouped with CYP94B1[47], an enzyme catalyzed the hydroxylation at C12 of jasmonyl-L-amino acid. Transcript/17284, transcript/17057, transcript/17636 and transcript/17854 were grouped in second clade. Transcript/17057 fell into a group with CYP82C4[48], an enzyme hydrolyzed xanthotoxin (8-methoxypsoralen) into 5-hydroxyxanthotoxin. Transcript/17284 grouped with CYP94B3[49] and it revealed that transcript/17284 may act as a hydroxylase. Transcript/17636 and transcript/17854 were divided into a group with uncharacterized CYP98A1.

According to CYP450 domain analysis, transcript/17636 and transcript/16971 contained the same CypX domain [50] with CYP81B1. P450-cyclo_AA_1 domain in transcript/16950 was also found within cytokinin trans-hydroxylase [51]. PLN00168 domain [52] was found in transcript/17824 and transcript/19854, and PLN02687 domain (a domain in flavonoid 3'-monooxygenase [53]) was also contained by transcript/16057, transcript/16577 and transcript/17057. Transcript/16777 and transcript/17001 contained PLN02774 domain, which was often found in brassinosteroid-6-oxidase [54]. Therefore, four DEGs (transcript/16057, transcript/ 16577 and transcript/17057 and transcript/17001) contained domains within CYP450 oxidases and the others contained domain within CYP450 hydroxylase

Among all ACT DEGs, transcript/18186 was grouped with ARE1, the enzyme encoding sterol O-acyltransferase [55]. Transcript/18214 fell into a subclade with LPAT2, an acyl-sn-glycerol-3-phosphate acyltransferase of Brassica oleracea [56]. Transcript/17792 was grouped with HCT2 which catalyzing the transfer of an acyl from p-coumaroyl-CoA to various acyl acceptors. Transcript/19132 and TSM1[57] were in a group which reveal that transcript/19132 was likely to be methyltransferase. Through domain analysis, transcript/17792 and transcript/18214 contained HXXXD domain that often found in BAHD family [58]. Transcript/19132 contained same domain PLN02177 with glycerol-3-phosphate acyltransferase [59]. Transcript/18186 and eukaryotic initiation factor 4B [60] had the same eIF-4B domain. Therefore, transcript/17792 and transcript/18214 were selected as candidates of ACT involved in azadirachtin A biosynthesis combined the prediction of phylogenetic analysis and domain prediction.

For all EST DEGs after phylogenetic analysis, transcript/19998 and transcript/16750 grouped with KAI2[61]. KAI2 was reported to be involved in seed germination and didn’t show activity as esterase. Transcript/19188 was divided into a subclade with A. thaliana GDL15, which belonged to GDSL [62]-like lipase/acylhydrolase superfamily and displayed hydrolytic activity with esters. Transcript/18100 was grouped with PME3, a pectinesterase catalyzed the hydrolysis of (1,4)-α-D-galacturonosyl methyl ester [63]. Transcript/19748 formed a tight subclade with TGL1, which was a sterol esterase mediating the hydrolysis of steryl esters [64]. Transcript/19882 and transcript/19697 were in a group with HIDH [65] and CXE18, respectively and these two enzymes shew activity to carboxylic esters. The conserved domain analysis of EST candidates presented transcript/19188 contained Ser-His-Asp (Glu) triad found in a SNGH plant lipase [66]. Transcript/19882, transcript/19697, transcript/19748 had AES domain [67] which contained by acetyl esterase/lipase. PAE domain [68] in transcript/18100 was also found in pectin acetylesterase while transcript/16750 and transcript/19998 contained the plant pectinesterase inhibitor domain PLN02201. Therefore, transcript/16750 and transcript/19998 were removed from the list of candidates because of the prediction from phylogenetic analysis and domain analysis.

Molecular Docking analysis of CYP450s

From the phylogenic analysis of CYP450s, five CYP450s were grouped with the members in CYP71 and CYP72 clades. Molecular docking of five CYP proteins (CYP16057 and CYP16577, CYP16950, CYP16777 and CYP17001) was performed with four triterpenoids (ligands) (Figure 5 and Figure S4). Analysis revealed that binding energy was lowest in case of CYP16057 docked with tirucalla-7,24-dien-3β-ol forming zero hydrogen bond. However, azadirone and nimbolide formed stable complexes with CYP16057 with one hydrogen bond with -7.20 and-6.40 kcal/mol of binding energies (Figure 5 and Table S7), respectively. The docking analysis for CYP16577 revealed that, among all the ligands, binding energy was lowest for tirucalla-7,24-dien-3β-ol and nimbin with -10.07 and -9.83 kcal/mol forming zero and two hydrogen bonds respectively. Docking of CYP16577 with triterpenoids showed interaction through only one hydrogen bond with binding energy of -9.42 and-9.16 kcal/mol for azadirone and nimbolide (Figure 5 and Table S8). The binding energy of CYP16777 docked with azadirone and tirucalla-7,24-dien-3β-ol was -9.96 and -9.75 kcal/mol forming one hydrogen bond respectively. However, nimbolide and nimbin formed stable complexes with three and two hydrogen bonds respectively, indicating that the conformation of nimbolide is best suitable for CYP16777 when number of hydrogen bonds formed between protein and ligand is set as criteria. The hydrogen bonds between nimbolide and CYP16777 are formed at PHE354, ARG355 and ARG419 with ligand moiety at different positions (Table S9) with varying bond length (Figure 5). The docking analysis for CYP16950 revealed that, among all the ligands, binding energy was lowest for tirucalla-7,24-dien-3β-ol and azadirone with -8.31 and -7.90 kcal/mol without forming hydrogen bond. However, nimbin and nimbolide formed stable complexes with CYP16950 with one hydrogen bond with -6.32 and-6.66 kcal/mol of binding energies (Figure 5 and Table S10). Docking of CYP17001 with triterpenoids showed interaction through only one hydrogen bond with the lowest binding energy of -10.11kcal/mol for tirucalla-7,24-dien-3β-ol. Azadirone and nimbolide respectively formed stable complexes with CYP17001 through two hydrogen bonds. The hydrogen bonds between nimbin and CYP17001 are formed at ARG355, ARG419 and GLY423 with the binding energy is -9.82kcal/mol (Figure 5 and Table S11).

Tirucalla-7,24-dien-3β-ol was confirmed to be the scaffold of azadirachtin A. Nimbin, nimbolide and azadirone are three important compounds isolated from neem. They were proposed as intermediates in azadirachtin A pathway. According to the docking of five CYP450s with them, CYP16057, CYP16577, CYP16950 and CYP17001 showed strongest binding with tirucalla-7,24-dien-3β-ol. CYP16777 was more easily to bind with azadirone. Three residues in CYP16777 and CYP17001 respectively formed stable hydrogen bonds with nimbolide and nimbin. All the docking results indicate the priority of reactions between five CYP450s and four ligands. It also provided a theoretical basis for further functional assay of these CYP450s. The residues in proteins forming hydrogen bond with ligands also provided sites for mutation when we need to improve catalytic ability of these CYP450s.

Measurement of unigenes expression in secondary metabolite pathways in neem

The expressions of unigenes involved in secondary metabolites including three terpenoids, two sterols and putative azadirachtin A downstream pathway were analyzed based on the KEGG annotation and FPKM (Figure 6). Among all unigenes, 13 unigenes were found to be related to MVA pathway and 38 unigenes related to MEP pathway (Table S5). Some of them ((unigenes encoding mevalonate kinase (MVK), 1-deoxy-D-xylulose-5-phosphate synthase (DXPS) and 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase (MDS)) were expressing higher in leaf and fruit. GPP (geranyl pyrophosphate), catalyzed by Geranyl diphosphate synthase (GPPS) was the common intermediate of monoterpenoids. Four enzymes involved in myrcene, limonene and terpineol biosynthetic pathways expressed high in flower and fruit. As for (E, E)-4, 8, 12-trimethyltrideca-1, 3, 7, 11-tetraene (TMTT) biosynthesis, unigenes encoding geranylgeranyl diphosphate synthase (GGPPS) and CYP82G1 respectively expressed higher in flower and stem.

Farnesyl diphosphate synthase (FDS) catalyzed formation of farnesyl pyrophosphate (FPP) from IPP and unigene encoding FDS expressed highest in fruit followed by in root and flower. Solavetivol pathway is one of sesquiterpenoid pathways detected and unigene encoding solavetivol synthase (CYP71D55) expressed higher in fruit and leaf. In triterpenoids biosynthesis, two FPPs formed to 2,3-oxidosqualene under continuously catalysis of squalene synthase (SQS) and squalene epoxidase (SQLE). SQS and SQLE encoded unigenes both expressed the highest in leaf. Unigenes encoding cycloartenol synthase (CAS1) was identified and its expression in different tissues was in the order of fruit> stem>root. Methylsterol monooxygenase (SMO1) and sterol-4alpha-carboxylate 3-dehydrogenase (NSDHL) expressed the highest in flower. Unigenes encoding delta14-sterol reductase (TM7SF2) and CYP51G1 were highly-expressed in leaf and fruit respectively.

The expression level of unigenes involved in putative azadirachtin A downstream pathway had been presented. Transcript/14449 encoding the first enzyme in azadirachtin A downstream pathway expressed highest in fruit followed by in leaf. After scaffold synthesis, the methyl group at C14 was removed. It was catalyzed by enzyme encoded by transcript/16198, which was highly-expressed in leaf. Alcohol at C3 was continuously oxidized into C3 ketone group by transcript/18725 and transcript/17679 and formed common compounds [69] isolated from Meliaceae family. These two unigenes expressed highest in root and fruit respectively. Next important step involved in azadirachtin A was the formation of furan ring, and two CYP450s catalyzing the formation were isolated from M. azedarach, and C. sinensis [18]. Two transcripts (transcript/16971 and transcript/16742) in our study were found to be homologs of the identified two CYP450s and might produce melianone [70] from the precursor. Transcript/16971 and transcript/16742 expressed highest in leaf and fruit, respectively. With some unknown enzymes, melianone was further modified into compound with furan ring and C7-OH. The insecticidal C7-hydroxylased compound [71] was esterized by transcript/18100 (higher-expressed in fruit and leaf) and further formed compounds like nimbin or nimbolide after some modifications. However, reactions between nimbin or nimbolide and azadirachtin A were still unclear.

After the bioinformatic analysis of these DEGs, two transcripts (transcript/18833 and transcript/19291) encoding ADH, twelve transcripts encoding CYP450, two ACT transcripts (transcript/17792 and transcript/18214) and five transcripts encoding EST were selected as candidates in azadirachtn A downstream pathway. Some DEGs were removed from the library after phylogenetic analysis and domain prediction, and some non-encoding DEGs were also deleted. DEGs were also selected although they were not higher expressed in fruit and leaf, for example transcript/16198, transcript/16742 and transcript/16950. Transcript/16198 was annotated to encode sterol 14-demethylase, which removes methyl group from C14 of sterol. Transcript/16742 encoded protein with 509 amino acids and it was 98% identical to MaCYP71BQ5(Figure S5). MaCYP71BQ5 is the CYP450 found to be involved in melianol formation. Researchers could only get the fragment of its homolog (AiCYP71BQ5) from neem while our transcript/16742 contained the complete ORF of AiCYP71BQ5.

As a kind of terpenoid, the increase of precursor leads to the higher production of terpenoid. In the case of artemisinic acid production, improvement of terpenoid precursor by engineering MVA pathway made 500 time increases in yield [72]. Thus, the up-regulated DEGs in MVA or MEP pathway and 2,3-oxidosqualene biosynthetic pathway would be used as building blocks in the construction of azadirachtin A precursor biosynthetic pathway in the future.

Although reaction types and key enzymes were partially proposed based on the structural differences between intermediates in our putative azadirachtin A pathway, some information was still missing. For instance, we couldn’t find the enzyme catalyzes the hydroxylation reaction at C7 site. What’s more, the order of reactions during the azadirachtin A downstream was not clear. Neither the numbers of reactions nor the catalysis type were characterized. These limitations of pathway lead to insufficient mining on neem transcriptome data. This might be one of the reasons for slow progress in azadirachtin pathway exploration although numerous available neem genome and transcriptome data.

In conclusion, the multi-tissue transcriptome analysis revealed five types of genes potentially involved in azadirachtin A downstream pathway and their respective transcription levels. It also indicated that the neem tree genome encodes a high number of terpene or limonoid biosynthetic genes. Finally, 22 unigenes encoding enzymes include OSC, ADH, CYP450, ACT and EST were selected as candidates involved in azadirachtin A downstream pathway. This is the first report on hybrid-seq transcriptome profiling analysis of A. indica. The obtained unigenes may provide a valid and diverse candidate pool for the study of the selective modification by the functional groups at triterpenoid or limonoids skeleton as well as the convergent evolution in secondary metabolism.

DEGs: differentially expressed genes; OSC: 2,3-oxidosqualene cyclase; ADH: alcohol dehydrogenase; CYP450: cytochrome P450; ACT: acyltransferase; EST: esterase; PacBio SMRT: Pacific Biosciences Single Molecule Real Time; FPKM: Fragments per kilobase per transcript per million mapped reads; Nr: non-redundant protein database, Nt: nucleotide sequences; COG: Clusters of Orthologous Groups of proteins; KEGG: Kyoto Encyclopedia of Genes and Genomes protein database; GO: Gene Ontology; qRT-PCR: quantitative real-time PCR; NGS: next-generation sequencing; IDI: Isopentenyl-diphosphate δ-isomerase; GGPS: Geranyl diphosphate synthase; FDS: Farnesyl diphosphate synthase; SQLE: Squalene epoxidase; G3P: 3-phosphate glyceraldehyde; MVA: mevalonate; MEP: methylerythritol phosphate; IPP: isopentenyl pyrophosphate; DMAPP: γ-dimethylallyl pyrophosphate; GPP: geranyl pyrophosphate; FPP: farnesyl pyrophosphate; DXP :1-Deoxy-D-xylulose 5-phosphate; MVP: 5-phosphomevalonate; MVPP: (R)-5-Diphosphomevalonate; GGPP: geranylgeranyl pyrophosphate; HMG-CoA: 3-hydroxy-3-methyl-glutaryl-CoA;PCME:2-phospho-4-2-C-methyl-D-erythritol; MECP: 2-C-methyl-D-erythritol 2,4-cyclodiphosphate; HMBDP: 1-hydroxy-2-methyl-2-butenyl 4-diphosphate; TMTT: (E,E)-4,8,12-trimethyltrideca-1,3,7,11-tetraene; DXPS: 1-deoxy-D-xylulose-5-phosphate synthase; MDS: 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase; HMGR: hydroxymethyl-glutaryl-CoA reductase; MVK: Mevalonate kinase; PMVK: Phosphomevalonate kinase; GPPS: Geranyl diphosphate synthase; MCS: Myrcene/ocimene synthase; LMS: (R)-limonene synthase; ATNS: (-)-alpha-terpineol synthase; CNS: 1,8-cineole synthase; TM7SF2: Delta14-sterol reductase; NSDHL: Sterol-4alpha-carboxylate 3-dehydrogenase; SMO1: Methylsterol monooxygenase 1; GGPPS: Geranylgeranyl diphosphate synthase; SQS: Squalene synthase; CAS1: Cycloartenol synthase.

Ethics approval and consent to participate

The locations of material collected here are neither privately owned lands nor protected areas. No specific permits were required for our research.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

We thank Zhiwei Wang and Huangying Shu (Hainan Key Laboratory for Sustainable Utilization of Tropical Bioresources, College of Horticulture, Hainan University) for neem identifying and sampling.

Funding

This work was supported by National Key R&D Program of China(2017YFD0201400) and the Fundamental Research Funds for the Central Universities. The funders did not have any role in the design of the study, collection, analysis, or interpretation of data or in the writing of the manuscript.

Availability of data and materials

All data generated or analyzed during the current study are included in this article and its supplementary information files. The raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database under BioProject ID PRJNA590058.

Author’s contributions

HW and YH conceived this study and designed the experiments. HW conducted all the experiments and analyzed the data. HW wrote the manuscript, NW and YH reviewed and edited the manuscript. All authors read and approved the final manuscript.

Hummel H, Langner S, Leithold G, Schmutterer H. NEEM: UNUSUALLY VERSATILE PLANT GENUS AZADIRACHTA WITH MANY USEFUL AND SO FAR INSUFFICIENTLY EXPLOITED PROPERTIES FOR AGRICULTURE, MEDICINE, AND INDUSTRY. Commun Agric Appl Biol Sci. 2014; 79:211-228.
Agbo B, And A, Ajaba M. A REVIEW ON THE USE OF NEEM (Azadirachta indica) AS A BIOPESTICIDE. Journal of Biopesticide and Environment. 2015; 2.
Oulhaci CM, Denis B, Kilani-Morakchi S, Sandoz J-C, Kaiser L, Joly D, et al. Azadirachtin effects on mating success, gametic abnormalities and progeny survival in Drosophila melanogaster (Diptera). Pest Manag Sci. 2018; 74(1):174-180.
P&S Market Research. Global Neem Extract Market Size, Share, Development, Growth and demand Forecast to 2022 [https://www.psmarketresearch.com/market-analysis/neem-extract-market]
Ambrosino P, Fresa R, Fogliano V, Monti SM, Ritieni A. Extraction of Azadirachtin A from Neem Seed Kernels by Supercritical Fluid and Its Evaluation by HPLC and LC/MS. J Agric Food Chem. 1999; 47(12):5252-5256.
Bilton JN, Broughton HB, Jones PS, Ley SV, Rzepa HS, Sheppard RN, et al. An x-ray crystallographic, mass spectroscopic, and NMR study of the limonoid insect antifeedant azadirachtin and related derivatives. Tetrahedron. 1987; 43(12):2805-2815.
Veitch GE, Beckmann E, Burke BJ, Boyer A, Maslen SL, Ley SV. Synthesis of Azadirachtin: A Long but Successful Journey. Angew Chem, Int Ed. 2007; 46(40):7629-7632.
Ekong DEU, Ibiyemi SA, Olagbemi EO. The meliacins (limonoids). Biosynthesis of nimbolide in the leaves of Azadirachta indica. J Chem Soc D. 1971(18):1117-1118.
Narnoliya LK, Rajakani R, Sangwan NS, Gupta V, Sangwan RS. Comparative transcripts profiling of fruit mesocarp and endocarp relevant to secondary metabolism by suppression subtractive hybridization in Azadirachta indica (neem). Mol Biol Rep. 2014; 41(5):3147-3162.
Rajakani R, Narnoliya L, Sangwan NS, Sangwan RS, Gupta V. Subtractive transcriptomes of fruit and leaf reveal differential representation of transcripts in Azadirachta indica. Tree Genetics & Genomes. 2014; 10(5):1331-1351.
Krishnan NM, Pattnaik S, Jain P, Gaur P, Choudhary R, Vaidyanathan S, et al. A draft of the genome and four transcriptomes of a medicinal and pesticidal angiosperm Azadirachta indica. BMC Genom. 2012; 13:464-464.
Krishnan NM, Jain P, Gupta S, Hariharan AK, Panda B. An Improved Genome Assembly of Azadirachta indica A. Juss. G3 (Bethesda). 2016; 6(7):1835-1840.
Wang S, Zhang H, Li X, Zhang J. Gene expression profiling analysis reveals a crucial gene regulating metabolism in adventitious roots of neem (Azadirachta indica). RSC Adv. 2016; 6(115):114889-114898.
Pandreka A, Dandekar DS, Haldar S, Uttara V, Vijayshree SG, Mulani FA, et al. Triterpenoid profiling and functional characterization of the initial genes involved in isoprenoid biosynthesis in neem (Azadirachta indica). BMC Plant Biol. 2015; 15:214-214.
Krishnan N, Pattnaik S, Sa D, K Hariharan A, Gaur P, Chaudhary R, et al. De novo sequencing and assembly of Azadirachta indica fruit transcriptome. Curr Sci. 2011; 101:1553.
Bhambhani S, Lakhwani D, Gupta P, Pandey A, Dhar YV, Kumar Bag S, et al. Transcriptome and metabolite analyses in Azadirachta indica: identification of genes involved in biosynthesis of bioactive triterpenoids. Sci Rep. 2017; 7(1):5043.
Wang Y, Chen X, Wang J, Xun H, Sun J, Tang F. Comparative analysis of the terpenoid biosynthesis pathway in Azadirachta indica and Melia azedarach by RNA-seq. Springerplus. 2016; 5(1):819-819.
Hodgson H, De La Peña R, Stephenson MJ, Thimmappa R, Vincent JL, Sattely ES, et al. Identification of key enzymes responsible for protolimonoid biosynthesis in plants: Opening the door to azadirachtin production. Proc Natl Acad Sci. 2019; 116(34):17096-17104.
Kurimoto S-i, Takaishi Y, Ahmed FA, Kashiwada Y. Triterpenoids from the fruits of Azadirachta indica (Meliaceae). Fitoterapia. 2014; 92:200-205.
Song L, Wang J, Gao Q, Ma X, Wang Y, Zhang Y, et al. Simultaneous determination of five azadirachtins in the seed and leaf extracts of Azadirachta indica by automated online solid-phase extraction coupled with LC–Q-TOF–MS. Chem Cent J. 2018; 12(1):85.
Sundaram KMS. Azadirachtin biopesticide: A review of studies conducted on its analytical chemistry, environmental behaviour and biological effects. J Environ Sci Health B. 1996; 31(4):913-948.
Rhoads A, Au KF. PacBio Sequencing and Its Applications. Genomics Proteomics Bioinformatics. 2015; 13(5):278-289.
Bashir A, Klammer A, Robins WP, Chin C-S, Webster D, Paxinos E, et al. A hybrid approach for the automated finishing of bacterial genomes. Nat Biotechnol. 2012; 30(7):701-707.
Koren S, Schatz MC, Walenz BP, Martin J, Howard JT, Ganapathy G, et al. Hybrid error correction and de novo assembly of single-molecule sequencing reads. Nat Biotechnol. 2012; 30(7):693-700.
Goodwin S, Gurtowski J, Ethe-Sayers S, Deshpande P, Schatz MC, McCombie WR. Oxford Nanopore Sequencing, Hybrid Error Correction, and de novo assembly of a Eukaryotic Genome. bioRxiv. 2015:013490.
Xu J, Wang X-y, Guo W-z. The cytochrome P450 superfamily: Key players in plant development and defense. J Integr Agric. 2015; 14(9):1673-1686.
Paddon CJ, Keasling JD. Semi-synthetic artemisinin: a model for the use of synthetic biology in pharmaceutical development. Nat Rev Microbiol. 2014; 12(5):355-367.
Legrand G, Delporte M, Khelifi C, Harant A, Vuylsteker C, Mörchen M, et al. Identification and Characterization of Five BAHD Acyltransferases Involved in Hydroxycinnamoyl Ester Metabolism in Chicory. Front Plant Sci. 2016; 7(741).
Au KF, Underwood JG, Lee L, Wong WH. Improving PacBio long read accuracy by short read alignment. PLoS One. 2012; 7(10):e46679-e46679.
Hu Z, Li G, Sun Y, Niu Y, Ma L, He B, et al. Gene transcription profiling of Aspergillus oryzae 3.042 treated with ergosterol biosynthesis inhibitors. Braz J Microbiol. 2019; 50(1):43-52.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987; 4(4):406-425.
Kelley LA, Mezulis S, Yates CM, Wass MN, Sternberg MJE. The Phyre2 web portal for protein modeling, prediction and analysis. Nat Protoc. 2015; 10(6):845-858.
Pandey V, Dhar YV, Gupta P, Bag SK, Atri N, Asif MH, et al. Comparative interactions of withanolides and sterols with two members of sterol glycosyltransferases from Withania somnifera. BMC Bioinform. 2015; 16(1):120-120.
Zhang W, Chen J, Keyhani NO, Zhang Z, Li S, Xia Y. Comparative transcriptomic analysis of immune responses of the migratory locust, Locusta migratoria, to challenge by the fungal insect pathogen, Metarhizium acridum. BMC Genom. 2015; 16:867-867.
Everaert C, Luypaert M, Maag JLV, Cheng QX, Dinger ME, Hellemans J, et al. Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Sci Rep. 2017; 7(1):1559-1559.
Pombo MA, Zheng Y, Fei Z, Martin GB, Rosli HG. Use of RNA-seq data to identify and validate RT-qPCR reference genes for studying the tomato-Pseudomonas pathosystem. In: Sci Rep. vol. 7; 2017: 44905.
Shiroguchi K, Jia TZ, Sims PA, Xie XS. Digital RNA sequencing minimizes sequence-dependent bias and amplification noise with optimized single-molecule barcodes. Proc Natl Acad Sci U S A. 2012; 109(4):1347-1352.
Ramsköld D, Luo S, Wang Y-C, Li R, Deng Q, Faridani OR, et al. Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells. Nat Biotechnol. 2012; 30(8):777-782.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE Guidelines: Minimum Information for Publication of Quantitative Real-Time PCR Experiments. Clin Chem. 2009; 55(4):611-622.
Kim S-J, Kim M-R, Bedgar DL, Moinuddin SGA, Cardenas CL, Davin LB, et al. Functional reclassification of the putative cinnamyl alcohol dehydrogenase multigene family in Arabidopsis. Proc Natl Acad Sci U S A. 2004; 101(6):1455-1460.
Achkor H, Díaz M, Fernández MR, Biosca JA, Parés X, Martínez MC. Enhanced formaldehyde detoxification by overexpression of glutathione-dependent formaldehyde dehydrogenase from Arabidopsis. Plant Physiol. 2003; 132(4):2248-2255.
Jin Y, Zhang C, Liu W, Qi H, Chen H, Cao S. The cinnamyl alcohol dehydrogenase gene family in melon (Cucumis melo L.): bioinformatic analysis and expression patterns. PLoS One. 2014; 9(7):e101730-e101730.
Hoh F, Pons JL, Gautier MF, de Lamotte F, Dumas C. Structure of a liganded type 2 non-specific lipid-transfer protein from wheat and the molecular basis of lipid binding. Acta Crystallogr D Biol Crystallogr. 2005; 61(Pt 4):397-406.
Kavanagh KL, Jörnvall H, Persson B, Oppermann U. Medium- and short-chain dehydrogenase/reductase gene and protein families : the SDR superfamily: functional and structural diversity within a family of metabolic and regulatory enzymes. Cell Mol Life Sci. 2008; 65(24):3895-3906.
Nagegowda DA, Gupta P. Advances in biosynthesis, regulation, and metabolic engineering of plant specialized terpenoids. Plant Sci. 2020; 294:110457.
Zheng X, Li P, Lu X. Research advances in cytochrome P450-catalysed pharmaceutical terpenoid biosynthesis in plants. J Exp Bot. 2019; 70(18):4619-4630.
Koo AJ, Thireault C, Zemelis S, Poudel AN, Zhang T, Kitaoka N, et al. Endoplasmic reticulum-associated inactivation of the hormone jasmonoyl-L-isoleucine by multiple members of the cytochrome P450 94 family in Arabidopsis. J Biol Chem. 2014; 289(43):29728-29738.
Rajniak J, Giehl RFH, Chang E, Murgia I, von Wirén N, Sattely ES. Biosynthesis of redox-active metabolites in response to iron deficiency in plants. Nat Chem Biol. 2018; 14(5):442-450.
Heitz T, Widemann E, Lugan R, Miesch L, Ullmann P, Desaubry L, et al. Cytochromes P450 CYP94C1 and CYP94B3 catalyze two successive oxidation steps of plant hormone Jasmonoyl-isoleucine for catabolic turnover. J Biol Chem. 2012; 287(9):6296-6306.
Cryle MJ, Bell SG, Schlichting I. Structural and Biochemical Characterization of the Cytochrome P450 CypX (CYP134A1) from Bacillus subtilis: A Cyclo-l-leucyl-l-leucyl Dipeptide Oxidase. Biochemistry. 2010; 49(34):7282-7296.
Takei K, Yamaya T, Sakakibara H. Arabidopsis CYP735A1 and CYP735A2 Encode Cytokinin Hydroxylases That Catalyze the Biosynthesis of trans-Zeatin. The Journal of biological chemistry. 2004; 279:41866-41872.
Lam PY, Liu H, Lo C. Completion of Tricin Biosynthesis Pathway in Rice: Cytochrome P450 75B4 Is a Unique Chrysoeriol 5′-Hydroxylase. Plant Physiol. 2015; 168(4):1527-1536.
Rice Annotation P, Itoh T, Tanaka T, Barrero RA, Yamasaki C, Fujii Y, et al. Curated genome annotation of Oryza sativa ssp. japonica and comparative genome analysis with Arabidopsis thaliana. Genome Res. 2007; 17(2):175-183.
Shimada Y, Fujioka S, Miyauchi N, Kushiro M, Takatsuto S, Nomura T, et al. Brassinosteroid-6-Oxidases from Arabidopsis and Tomato Catalyze Multiple C-6 Oxidations in Brassinosteroid Biosynthesis. Plant Physiol. 2001; 126(2):770-779.
Ivashov VA, Zellnig G, Grillitsch K, Daum G. Identification of triacylglycerol and steryl ester synthases of the methylotrophic yeast Pichia pastoris. Biochim Biophys Acta. 2013; 1831(6):1158-1166.
Kim HU, Li Y, Huang AHC. Ubiquitous and Endoplasmic Reticulum–Located Lysophosphatidyl Acyltransferase, LPAT2, Is Essential for Female but Not Male Gametophyte Development in Arabidopsis. The Plant Cell. 2005; 17(4):1073-1089.
Fellenberg C, Milkowski C, Hause B, Lange PR, Vogt T. Tapetum-specific location of a cation-dependent O-methyltransferase in Arabidopsis thaliana. Plant Journal. 2008; 56(1):132-145.
Panikashvili D, Shi JX, Schreiber L, Aharoni A. The Arabidopsis DCR encoding a soluble BAHD acyltransferase is required for cutin polyester formation and seed hydration properties. Plant Physiol. 2009; 151(4):1773-1789.
Yu J, Loh K, Song Z-y, Yang H-q, Zhang Y, Lin S. Update on glycerol-3-phosphate acyltransferases: the roles in the development of insulin resistance. Nutrition & Diabetes. 2018; 8(1):34.
Metz AM, Wong KCH, Malmström SA, Browning KS. Eukaryotic Initiation Factor 4B from Wheat and Arabidopsis thaliana Is a Member of a Multigene Family. Biochem Biophys Res Commun. 1999; 266(2):314-321.
Guo Y, Zheng Z, La Clair JJ, Chory J, Noel JP. Smoke-derived karrikin perception by the α/β-hydrolase KAI2 from Arabidopsis. Proc Natl Acad Sci U S A. 2013; 110(20):8284-8289.
Chepyshko H, Lai C-P, Huang L-M, Liu J-H, Shaw J-F. Multifunctionality and diversity of GDSL esterase/lipase gene family in rice (Oryza sativa L. japonica) genome: new insights from bioinformatics analysis. BMC Genom. 2012; 13:309-309.
Christensen TMIE, Nielsen JE, Kreiberg JD, Rasmussen P, Mikkelsen JD. Pectin methyl esterase from orange fruit: characterization and localization by in-situ hybridization and immunohistochemistry. Planta. 1998; 206(4):493-503.
Köffel R, Tiwari R, Falquet L, Schneiter R. The Saccharomyces cerevisiae YLL012/YEH1, YLR020/YEH2, and TGL1 genes encode a novel family of membrane-anchored lipases that are required for steryl ester hydrolysis. Mol Cell Biol. 2005; 25(5):1655-1668.
Akashi T, Aoki T, Ayabe S-I. Molecular and biochemical characterization of 2-hydroxyisoflavanone dehydratase. Involvement of carboxylesterase-like proteins in leguminous isoflavone biosynthesis. Plant Physiol. 2005; 137(3):882-891.
Mølgaard A, Kauppinen S, Larsen S. Rhamnogalacturonan acetylesterase elucidates the structure and function of a new family of hydrolases. Structure. 2000; 8(4):373-383.
Pereira EO, Tsang A, McAllister TA, Menassa R. The production and characterization of a new active lipase from Acremonium alcalophilum using a plant bioreactor. Biotechnol Biofuels. 2013; 6:111-111.
Philippe F, Pelloux J, Rayon C. Plant pectin acetylesterase structure and function: new insights from bioinformatic analysis. BMC Genom. 2017; 18(1):456.
Paal C. Ueber die Derivate des Acetophenonacetessigesters und des Acetonylacetessigesters. Ber Dtsch Chem Ges. 1884; 17(2):2756-2767.
Polonsky J, Varon Z, Rabanal RM, Jacquemin H. 21,20-Anhydromelianone and Melianone from Simarouba amara (Simaroubaceae); Carbon-13 NMR Spectral Analysis of Δ7-Tirucallol-Type Triterpenes. Isr J Chem. 1977; 16(1):16-19.
Huidana F, Conghaia Z, Shengjiaoa Y, Jun L. Advances of synthesis and structure modification and bioactivity of azadirachtin. Chinese J Org Chem. 2009; 29:20-33.
Paddon CJ, Westfall PJ, Pitera DJ, Benjamin K, Fisher K, McPhee D, et al. High-level semi-synthetic production of the potent antimalarial artemisinin. Nature. 2013; 496:528.

Table S1. Concentration of metabolites in leaf, bark and seed extract of A. indica.

Table S2. Summary of functional annotations of A. indica.

Table S3. KEGG annotation of unigenes from A. indica.

Table S4. Pathways and number of unigenes related to secondary metabolites in A. india.

Table S5. Discovery and expression of unigenes involved in terpenoid backbone biosynthesis in A. indica.

Table S6. Unigenes involved in the sesquiterpenoid and triterpenoid biosynthesis in A. indica.

Table S7. Details of molecular docking of CYP16057 protein with ligands.

Table S8. Details of molecular docking of CYP16577 protein with ligands.

Table S9. Details of molecular docking of CYP16777 protein with ligands.

Table S10. Details of molecular docking of CYP16950 protein with ligands.

Table S11. Details of molecular docking of CYP17001 protein with ligands.

Additional file 1: Details of proteins used for phylogenetic analysis of OSC, ADH, CYP450, ACT and EST candidates.

Additional file 2: Genes primers used for qRT-PCR analysis and full-length cloning PCR.

Figure S1. The work-flow of gene mining. Letters L, F, R, S and f represent leaf, fruit, root, stem and flower, respectively.

Figure S2. Quantitative Real-Time PCR (qRT-PCR) validation of selected unigenes from fruit, leaf, root, stem and flower tissues of A. indica. The relative expression level of each selected gene was determined by 2-^𝛥𝛥CT. Experiments have been conducted in triplicates.

Figure S3. The alignment result of transcript/14449 and AiOSC1. Sequence marked with “*” represented transcript/14449 from our data.

Figure S4. Molecular docking analysis of CYP450s putatively involved in azadirachtin A biosynthesis. Interaction of five complete protein (CYP16057, CYP16577, CYP16777, CYP16950 and CYP7001) with four ligands (tirucalla-7,24-dien-3β-ol, azadirone, nimbin and nimbolide) are shown in the figure. Five CYP450s and four ligands are represented in different colors. Complete block represents the whole protein ligand complex.

Figure S5. The amino acid sequence alignment of transcript/16742, AiCYP71BQ5 and MaCYP71BQ5. Same sequences were shown in white with red background.

Table 2. Summary of Illumina Hiseq and PacBio data output quality and assembled sequences of six libraries of Azadirachtin indica

Sequencing Platform	Sample	Reads Number (M)	Clean Bases (G)	Q20 (%)	Q30 (%)	Total number of unigene	Total length of unigene (bp)	Mean length of unigene (bp)	N50	GC (%)
Ilumina HiSeq	Leaf	41.14	6.17	96.95	92.30	50394	82508501	1380	2201	40.64
	Flower	41.35	6.20	97.29	93.06	62426	95564612	1530	2335	40.13
	Stem	40.60	6.09	97.29	93.11	65762	98612970	1499	2372	40.39
	Fruit	40.59	6.09	97.45	93.39	66668	81799852	1226	2077	41.93
	Root	40.64	6.10	97.44	93.43	45459	57083335	1255	2100	41.10
Hiseq summary	—	—	—	—	—	113008	175268545	1550	2599	40.99
PacBio SMRT	Mixed tissue	6.75	11.28	—	—	22884	82035635	3584	5068	43.55
Calibrated PacBio by Illumina		—	—	—	—	20201	72872459	3607	5076	43.55

Q20 and Q30 on Illumina platform correspond to the predicted base call error rate of 1 % and 0.1%, respectively.

N50: The minimum contig length needed to cover 50% of the transcriptome.

Download PDF

Journal Publication

published 28 Oct, 2020

Read the published version in BMC Genomics →

Review #1 received at journal
18 Jun, 2020
Editorial decision: Major revision
18 Jun, 2020
Reviewer #2 agreed at journal
31 May, 2020
Review #2 received at journal
31 May, 2020
Reviewer #1 agreed at journal
25 May, 2020
Editor assigned by journal
24 May, 2020
Reviewers invited by journal
24 May, 2020
Submission checks completed at journal
23 May, 2020
Editor invited by journal
23 May, 2020

You are reading this older preprint version

Read the latest preprint version →

Multi-tissue transcriptome analysis using hybrid-seq revealed potential genes and biological pathways associated with azadirachtin A biosynthesis in neem (Azadirachtin indica)

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Methods

Results And Discussions

Conclusions

Abbreviations

Declarations

References

Additional Files

Table 2

Supplementary Files

Status:

Journal Publication

Version 2