Computational identification of sweet wormwood (Artemisia annua) microRNA and their mRNA targets.

Despite its efficacy against malaria, the relatively low yield (0.01%-0.8%) of artemisinin in Artemisia annua is a serious limitation to the commercialization of the drug. A better understanding of the biosynthetic pathway of artemisinin and its regulation by both exogenous and endogenous factors is essential to improve artemisinin yield. Increasing evidence has shown that microRNAs (miRNAs) play multiple roles in various biological processes. In this study, we used previously known miRNAs from Arabidopsis and rice against expressed sequence tag (EST) database of A. annua to search for potential miRNAs and their targets in A. annua. A total of six potential miRNAs were predicted, which belong to the miR414 and miR1310 families. Furthermore, eight potential target genes were identified in this species. Among them, seven genes encode proteins that play important roles in artemisinin biosynthesis, including HMG-CoA reductase (HMGR), amorpha-4,11-diene synthase (ADS), farnesyl pyrophosphate synthase (FPS) and cytochrome P450. In addition, a gene coding for putative AINTEGUMENTA, which is involved in signal transduction and development, was also predicted as one of the targets. This is the first in silico study to indicate that miRNAs target genes encoding enzymes involved in artemisinin biosynthesis, which may help to understand the miRNA-mediated regulation of artemisinin biosynthesis in A. annua.

Introduction

Malaria is a global health problem with more than 1 billion people living in areas at a high risk of the disease. Artemisinin combination therapies (ACTs) are the recommended treatment regimen for uncomplicated malaria caused by the Plasmodium falciparum parasite (. It has long been recognized that the problem of artemisinin resistance is best addressed by increasing access to ACTs (. This approach receives strong support from the global health community. However, there is growing concern that the supply chain will be unable to consistently produce high-quality artemisinin in the quantities that will be required (. ACT supply remains reliant on the agricultural production of artemisinin, which is a sesquiterpenoid synthesized in the glandular trichomes of the Chinese medicinal plant sweet wormwood (Artemisia annua L.) (. However, artemisinin concentration in A. annua is low, (in the range of 0.01%-0.8% per dry weight of tissue), which seriously limits the commercialization of the drug (. Therefore, improved varieties of A. annua for farmers in the developing countries would bring immediate benefits to the existing artemisinin supply chain by reducing production costs, stabilizing supplies, and improving growers’ confidence in the crop (.

Conventional breeding and genetic engineering approaches would allow construction of A. annua genotypes rich in artemisinin (. A better understanding of the biochemical pathway leading to the artemisinin synthesis and regulation by both exogenous and endogenous factors is essential for facilitating yield increase (. With the elucidation of the artemisinin biosynthetic pathway and identification of amorpha-4,11-diene synthase (ADS), which catalyses the first biosynthetic step in artemisinin biosynthesis, it is possible to explore unconventional alternate strategies that are economically viable for the commercial production of artemisinin. Two approaches appear promising. The first approach is to synthesize artemisinin from its simple precursor such as artemisinic acid via semi-synthetic route. Recently, scientific effort is being directed to develop a biological method to supply sufficient and reliable quantities of artemisinic acid, a direct precursor of artemisinin. For example, there was report on engineering of Saccharomyces cerevisiae to produce high titres (up to 100 mg/L) of artemisinic acid using an engineered mevalonate pathway, ADS, and a novel cytochrome P450 monooxygenase (CYP71AV1) from A. annua (. The second approach is to regulate key enzymes leading to increased artemisinin biosynthesis with metabolic engineering. Certainly, the latter strategy has provided some exciting results and further efforts may accelerate commercialization of this crucial drug. Zhang et al. showed that down-regulation of squalene synthase (SQS), a key enzyme of sterol pathway that is competitive with artemisinin biosynthetic pathway, by hairpin-RNA-mediated gene silencing in A. annua resulted in a 3-fold increase in artemisinin production (. Furthermore, study on the expression of genes involved in the terpene metabolism also indicated that SQS may significantly compete for farnesyl diphosphate (FDP) in artemisinin-producing tissues of A. annua (. In addition, higher artemisinin content was reported in induced tetraploid A. annua, which may result from the upregulated expression of some key enzyme genes related to artemisinin biosynthesis including ADS, farnesyl diphosphate synthase (FPS), HMG-CoA reductase (HMGR) and artemisinin metabolite-specific aldehyde hydrogenase 1 (. Therefore, a better understanding of the molecular mechanisms involved in the artemisinin biosynthesis and regulation will provide better strategies to develop new varieties with a higher content of artemisinin.

miRNAs are a large family of endogenous small RNAs containing ~22 nucleotides, which are derived from large precursors that are transcribed from non-protein-coding genes (. Plant miRNAs generally interact with their targets through perfect or near-perfect complementarity and repress translation 11, 12 or cleave targeted mRNAs (, thus negatively regulate the expression of their target genes (. Plant miRNAs target a large number of genes with functions in a range of development processes, including meristem cell identity (, leaf organ morphogenesis 11, 12, polarity and floral differentiation and development (. miRNAs are also reported to be involved in plant responses to biotic and environmental stresses (. The metabolite biosynthesis is also regulated by miRNAs (. A large number of miRNAs have so far been identified in various plant species. However, no miRNA from Asteraceae has been reported yet, reflecting a disparity between the important values of this plant family and insufficient molecular and genetic studies, including small RNA mediated gene regulation, in Asteraceae. To gain insight into miRNAs and their important regulatory functions in artemisinin biosynthetic pathway, we studied miRNA and their targets in A.annua using computational approach.

Computational or bioinformatics approach is one of the many approaches available for miRNA prediction (, which can discover miRNAs not only from species with full genomic and sufficient EST database available, but also from those with incomplete genomic information while with sufficient EST sequences available (. In addition, this approach is very useful for predicting miRNAs that usually cannot be detected by the direct cloning, particularly the low-abundance miRNAs. The computational approaches are based on homology search, gene search, neighbor stem-loop search, comparative genomic algorithm or phylogenetic shadowing (. Homology search, which can be further classified as genome-based search or EST-based search (, is based on conserved sequences and secondary structures and identifies miRNA genes by searching nucleotide databases using BLAST. Using homology search, orthologues of known miRNAs were revealed in different species, supporting that miRNAs are conserved in different species (. In addition, hundreds of new miRNAs were also identified using this method from the genomes of model species, such as Arabidopsis ( and rice (Oryza. sativa) (.

With the development of computational methods, several computer software programs have been developed to help identify plant potential miRNA target genes in mRNA sequences. Because almost all miRNAs show perfect or near-perfect complementarity with their targets in plants, it is much easier to predict miRNA targets using a BLAST search of mRNA database. More and more studies have shown the success of this powerful approach to select potential miRNA targets in mRNA sequences for experimental validation (.

Results and Discussion

Identification of potential miRNAs in A. annua

Most mature miRNAs are evolutionarily conserved from species to species within the plant kingdom, which facilitates the prediction of the existence of new miRNA orthologs or homologs in other plant species. In this study, we applied the comprehensive strategy to identify potential miRNAs in A. annua by searching EST against known miRNAs of a dicotyledonous plant Arabidopsis and a monocotyledonous rice.

Following the procedure depicted in Figure 1, 94,724 ESTs from A. annua were searched against 584 mature sequences of miRNAs from Arabidopsis and rice after removal of redundant sequences. In total, six potential miRNAs were predicted from A. annua (Figure 2). The six identified A. annua candidate miRNAs belong to two miRNA families. miR1310 family has one miRNA. miR414 family has five homologs with two from A. thaliana and three from O. sativa, respectively (Table 1).

Figure 1

Schematic diagram for searching potential miRNA genes in A. annua by identifying homologs of previously known plant miRNAs.

Figure 2

Mature and precursor sequences and the predicted stem and loop structures of newly identified miRNAs in A. annua. The mature miRNAs are highlighted in gray. A. miR414 (EY064998; homolog of A. thaliana); B. miR414 (EY057163; homolog of O. sativa); C. miR414 (EY082442; homolog of O. sativa); D. miR414 (EY107691; homolog of A. thaliana); E. miR1310 (EY057327; homolog of O. sativa); F. miR414 (EY083764; homolog of O. sativa).

Table 1

Newly identified miRNAs from ESTs of A. annua

miRNA	Reference Species	Gene ID	EST length (nt)	NM (nt)	LM (nt)	LP (nt)	Side	A + U (%)	MFE	MEFI
miR414	A. thaliana	EY107691	697	5	21	80	3′	55	14.7	0.41
	A. thaliana	EY064998	795	4	21	150	3′	52	46.1	0.64
miR414	O. sativa	EY057163	759	5	21	100	5′	70	26.2	0.87
	O. sativa	EY083764	678	2	21	90	3′	52.2	28.5	0.66
	O. sativa	EY082442	612	3	21	80	5′	48.8	22.5	0.55
miR1310	O. sativa	EY057327	381	3	23	80	5′	42.5	29.4	0.64

Note: NM, number of mismatch; LM, length of mature miRNAs; LP, length of precursor; MFE, minimal folding free energy; MFEI, minimal folding free energy index.

It is estimated that in plants, approximate 10,000 ESTs contain 1 miRNA. Therefore, the total of 94,724 ESTs in A. annua examined in this study may contain 9-10 miRNAs. However, in this study, only six miRNAs were predicted even with maximal 5-mismatches were allowed. The average length of ESTs is 654 nt and the longest is 795 nt, while most of miRNA precursors have 80-150 nt as identified by MirEval software (, suggesting that the EST may contain other element sequences in addition to miRNA precursor sequence (. The length of miRNA precursors in A. annua varied from 80 to 150 nt, with an average of 97 nt. The different sizes of the identified miRNAs within different families suggest that they may offer unique functions for regulation of miRNA biogenesis or gene expression (. The diversity of the identified miRNAs could be also found in the location of mature miRNA sequences. It is shown that the sequences of miR1310 and two members of miR414 from O. sativa family were located at the 5’ end of the miRNA precursors, while the other miR414 members were found at the 3’ end. Minimal folding free energy (MFE) is an important characteristic that determines the secondary structure of nucleic acids (DNA and RNA). The lower the MFE is, the higher the thermodynamically stable secondary structure of the corresponding sequence is (. The MFE index (MFEI) for each sequence was calculated as previously reported (. In this study, the MEFI values ranged from 0.41 to 0.87. miRNA precursor sequences have significantly higher MFEI value than other non-coding or coding RNAs. To avoid false calling of other RNAs as miRNA candidates, MFEI was also considered when predicting secondary structures (.

Identifying miRNAs using EST analysis ( has some advantages over other methods (. It has been suggested that most of the miRNAs predicted by EST analysis can be recovered by high-throughput deep sequencing (.

Although we have computationally identified six miRNAs, the number of miRNAs discovered is relatively small. A. annua belongs to Asteraceae, which is the largest plant family of vascular plants on Earth; comprising more than 23,000 genetically diverse and ecologically successful species (http://www.mobot.org/MOBOT/research/APweb/, Angiosperm Phylogeny Website, Version 9, 2011). Unfortunately, not a single miRNA from Asteraceae family has been deposited in the MiRbase (. We expect that as more miRNAs of this family are publicly available, more miRNAs will be identified in A. annua. Interestingly, a recently published study reported that 151 potentially conserved miRNAs belonging to 26 miRNA families were successfully identified and characterized using qPCR in 11 genus of Asteraceae (. In addition, another recent article reported 11 highly conserved miRNA precursors from 9 families using gene-oriented clusters of transcript sequences of A. annua with 88,174 UniGenes using a modified computational approach (. However, no miRNA targets were found among genes encoding enzymes involved in artemisinin biosynthesis, maybe because the conserved miRNAs preform evolutionarily stable functions, or miRNA-mediated regulation of artemisinin synthesis could be exerted primarily by novel or clade-specific miRNAs.

Prediction of potential targets of putative miRNAs in A. annua

Gaining insight into the miRNA targets will help us understand the spectrum of miRNA regulation and elucidate the functional importance of miRNAs. miRNAs may directly target transcription factors which affect plant development and specific genes which control metabolism as well (. To identify potential regulatory targets, we first searched mRNA database in A. annua and screened for mRNAs complementary to the six miRNAs with less than 4 mismatches. Gaps, G-U and other non-canonical pairs were not allowed and considered as mismatches. By screening against mRNA sequences of A. annua using the six newly identified miRNAs, we found 8 target genes complementary with less than 4 mismatches. Interestingly, one miRNA can be complementary to more than one regulatory target (Table 2). For example, six sequences were detected as targets of miR414 of A. annua.

Table 2

List of the potential targets of newly identified miRNAs in A. annua

miRNA	Gene ID	Target Gene ID	Target Protein	Target Function
miR414	EY107691	GQ468547	Putative AINTEGUMENTA	Transcription regulation
	EY064998	DQ363131	Putative flavonoid 3′-hydroxylase cytochrome P450	Oxidation-reduction process
miR1310	EY057327	AJ001539	EPS	Lyase activity
miR414	EY057163	CAB56499	Putative SQC	Lyase activity
	EY083764	FJ432667, AF327527	ADS	Biosynthesis of artemisinin
	EY082442	GQ420346	FPS	Biosynthesis of cholesterol, isoprene, lipid, steroid and sterol
		U14625	HMGR	Oxidation-reduction process
		DQ363131	Putative flavonoid 3′-hydroxylase cytochrome P450	Oxidation-reduction process

Note: EPS, epi-cedrol synthase; SQC, sesquiterpene cyclase; ADS, amorpha-4,11-diene synthase; FPS: farnesyl pyrophosphate synthase; HMGR, HMG-CoA reductase.

One target gene identified for rice miR414 homolog in A. annua (EY082442) is FPS (Table 2). FPS catalyzes two consecutive condensation reactions to produce FDP, which is the starting point of a large variety of essential isoprenoid end products, including artemisinin 39, 40. It is clear that the first dedicated step in the biosynthesis of artemisinin is the cyclization of FDP to form amorpha-4,11-diene catalyzed by ADS, one of the sesquiterpene cyclases (SQCs) (. The cyclase reaction establishes an important stereochemical framework upon which all other chemical modifications take place (. Interestingly, another rice miR414 homolog in A. annua, EY083764, is predicted to target ADS (Table 2). Modification of the amorpha-4,11-diene carbon skeleton to produce artemisinin acid was thought to involve a cytochrome P450 enzyme leading to the production of artemisinic alcohol, which could then be oxidized twice by either cytochrome P450 enzymes or dehydrogenases to yield artemisinic acid (. In our study, we found that both EY082442 and another Arabidopsis miR414 homolog in A. annua (EY064998) were predicted to target putative flavonoid 3′-hydroxylase cytochrome P450 (DQ363131).

In addition, EY082442 was also predicted to target HMGR (U14625), which is a rate-limiting enzyme of the mevalonate pathway. Recently it has been shown that HMGR expression limits artemisinin formation in A. annua (. Overexpression of ADS and HMGR led to significant increase in the artemisinin yield from the transgenic A. annua (.

A third rice miR414 homolog of A. annua, EY 57163, targets a putative SQC (CAB56499), while epi-cedrol synthase (EPS, AJ001539), an identified SQC converting FDP to 8-epicedrol, is the only target of miR1310 homolog of A. annua, EY057327 (Table 2). A diagram depicting the involvement of miRNAs and their targets was shown in Figure S1.

Other than the enzymes involved in production of secondary metabolites, our computational result demonstrates that A. annua miR414 (EY107691) targets the mRNA encoding putative AINTEGUMENTA protein (GQ468547) (Table 2), which is a transcription activator that recognizes and binds to the DNA consensus sequence 5′-CAC[AG]N[AT]TNCCNAN G-3′ (. AINTEGUMENTA is required for the initiation and growth of ovules integumenta, and also involved in organ initiation and development, including floral organs (. Interestingly, a positive correlation was noticed between the plant age and artemisinin yield (. Furthermore, artemisinin is present in high concentration in either flowers or leaves but low or zero in stems and roots (. For A. annua, the highest artemisinin concentration has been reported in leaves and flowers during full bloom stage, in comparison to the pre- and post-flowering stages 49, 50.

Conclusion

A. annua has received increasing attention due to its ability to produce artemisinin, which today is widely used for treatment of malaria. In addition to its anti-malarial properties, artemisinin is cytotoxic for cancer cells. Recent reports demonstrate that artemisinin inhibits the secretion and gene expression of tumor necrosis factor (TNF)-α, interleukin (IL)-1β, and IL-6 in a dose-dependent manner (. The unfortunately low yield of artemisinin and a worldwide shortage of the drug has led to intense research in order to increase the yield of this sesquiterpenoid. Our study is the first in silico study to identify miRNAs and their targets in A. annua, which we hope could help to better understand miRNA-mediated regulation of genes related to artemisinin biosynthesis.

In summary, in the present study, we predicted six miRNAs conserved in A. annua. Furthermore, eight potential target genes were predicted, with functions in a variety of biological processes, including artemisinin biosynthesis, signal transduction and development. Most of the targets are unique to A. annua genome, which encode the enzymes associated with the artemisinin biosynthesis pathway. Interestingly, we also identify one target coding for putative AINTEGUMENTA, a transcription factor involved in developmental process, especially for floral organs, which is coincident with the higher artemisinin content during flowering. The result from the computational prediction will be useful to guide experimental design for biological verification. The next major steps, therefore, are to experimentally analyze the functional categories suggested by our computational approach, determine the analogous molecular functions amongst divergent plant species, and further elucidate any significant correlations between the miRNAs and their target genes. Hopefully all these efforts would help make more artemisinin available at lower costs for more people in the Third World, so people who suffer most from malaria can benefit more from this valuable and effective drug.

Materials and Methods

Databases of miRNAs, ESTs, and mRNA sequences

To search potential miRNAs, a total of previously known 674 miRNAs and their precursor sequences from A. thaliana and O. sativa were obtained from miRNA Registry Database (Release 16.0, October 2010; http://www.mirbase.org/) (. These miRNAs were defined as the reference set of miRNA sequences. We have referred to the previous work on computational prediction of miRNAs by Zhang et al. (. To avoid the redundant or overlapping miRNAs, the repeated sequences of miRNAs within the above species were removed and the remaining 584 sequences were used as query sequences for BLAST search. A. annua EST and mRNA databases were obtained from the National Center for Biotechnology Information (NCBI) GenBank nucleotide databases (http://ftp.ncbi.nlm.nih.gov).

Availability of software

Comparative software BLAST-2.2.14 was used from NCBI GenBank. MFOLD 3.1 from website (http://www.bioinfo.rpi.edu/applications/mfold/rna/form1.cgi) was used online to analyze secondary structure of RNAs. MirEval (http://tagc.univ-mrs.fr/mireval) was used to predict miRNA precursors (. BLASTx from NCBI (http://www.ncbi.nlm.nih.gov) was used to analyze potential targets of miRNAs.

Prediction of miRNAs

Procedure for searching potential miRNAs in A. annua is shown in Figure 2. We used the method described by Zhang et al. ( with some modifications. Briefly, the previously known miRNAs in A. thaliana and O. sativa were screened out, and the redundant sequences were removed. The remaining miRNA sequences were subjected to BLAST search for A. annua miRNA homologs against EST databases.

The mature sequences of all miRNAs from Arabidopsis and rice were subjected to BLASTn search in the A. annua EST databases using BLASTn 2.2.9. The adjusted BLASTn parameter settings were as follows: expect values were set at 1,000; low complexity was chosen as the sequence filter; the number of descriptions and alignments was raised to 1,000. The default word-match size between the query and database sequences was 7. RNA sequences were considered as miRNA candidates only if they fit the following criteria: (1) at least 18 nt length were adopted between the predicted mature miRNAs and (2) allowed to have 0-5 nt mismatches in sequence with all previously known plant mature miRNAs (. The ESTs that closely match the previously known plant mature miRNAs were included in the set of miRNA candidates and used for additional characterization based on the following criteria: (1) the entire EST sequence was selected to predict the secondary structures and to screen for miRNA precursor sequences; (2) the selected ESTs were further compared with each other to eliminate redundancies; and (3) these sequences were subjected to evaluation for miRNA precursor prediction properties using mirEval software (. These precursor sequences were used for BLASTx analysis for removing the protein-coding sequences and retaining only the non-protein-coding sequences.

Prediction of secondary structure

Precursor sequences of these potential miRNA homologs were subjected to hairpin structure predictions using the Zuker folding algorithm with Mfold-3.1. The following parameters were used in predicting the secondary structures: (1) linear RNA sequence; (2) folding temperatures fixed at 37°C; ionic conditions of 1 M NaCl without divalent ions; (3) percent suboptimality number of 5; (4) maximum interior/bulge loop size of 30; (5) the grid lines in energy dot plot turned on. All other parameters were set with default values. In brief, the following criteria were applied in designating the RNA sequence as an miRNA homolog as described by Wang et al. (: (1) the sequence could fold into an appropriate stem-loop hairpin secondary structure; (2) the small RNA sits in one arm of the hairpin structure; (3) no more than 6 mismatches are between the predicted mature miRNA sequence and its opposite miRNA (miRNA*) sequence in the secondary structure; (4) no loop or break is in the miRNA or miRNA* sequences, and (5) predicted secondary structure has higher MFEI and negative MFE.

The MFEI was calculated using the following equation (:

MEFI= [(MEF/length of the RNA sequence) × 100] / (G+C)%

MFE denotes the negative folding free energies (ΔG).

Prediction of mRNA targets of miRNAs

Previous study has shown that most known plant miRNAs bind to the protein-coding region of their mRNA targets with perfect or nearly-perfect sequence complementarity, and degrade the target mRNA in a way similar to RNA interference 10, 12. This suggests a powerful approach to predict miRNA targets in plants by simply using homology search (Figure 3). In this study, we used homology search to predict miRNA targets in A. annua. The number of allowed mismatches at complementary sites between miRNA sequences and potential mRNA targets was no more than 4 and no gaps were allowed at the complementary sites.

Figure 3

Schematic diagram for searching potential target mRNA of miRNAs by blasting mRNA database of A. annua with newly identified miRNA sequences.

Authors’ contributions

RKM conceived the project. AP and RKM collected the data and conducted the computational analysis. NB and PKN supervised the work. AP and RKM interpreted the data and RKM prepared the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors have declared that no competing interests exist.