key: cord-308764-9z4qcoqz
authors: Wei, Lin; Li, Shenghua; Liu, Shenggui; He, Anna; Wang, Dan; Wang, Jie; Tang, Yulian; Wu, Xianjin
title: Transcriptome Analysis of Houttuynia cordata Thunb. by Illumina Paired-End RNA Sequencing and SSR Marker Discovery
date: 2014-01-02
journal: PLoS One
DOI: 10.1371/journal.pone.0084105
sha: 
doc_id: 308764
cord_uid: 9z4qcoqz

BACKGROUND: Houttuynia cordata Thunb. is an important traditional medical herb in China and other Asian countries, with high medicinal and economic value. However, a lack of available genomic information has become a limitation for research on this species. Thus, we carried out high-throughput transcriptomic sequencing of H. cordata to generate an enormous transcriptome sequence dataset for gene discovery and molecular marker development. PRINCIPAL FINDINGS: Illumina paired-end sequencing technology produced over 56 million sequencing reads from H. cordata mRNA. Subsequent de novo assembly yielded 63,954 unigenes, 39,982 (62.52%) and 26,122 (40.84%) of which had significant similarity to proteins in the NCBI nonredundant protein and Swiss-Prot databases (E-value <10(−5)), respectively. Of these annotated unigenes, 30,131 and 15,363 unigenes were assigned to gene ontology categories and clusters of orthologous groups, respectively. In addition, 24,434 (38.21%) unigenes were mapped onto 128 pathways using the KEGG pathway database and 17,964 (44.93%) unigenes showed homology to Vitis vinifera (Vitaceae) genes in BLASTx analysis. Furthermore, 4,800 cDNA SSRs were identified as potential molecular markers. Fifty primer pairs were randomly selected to detect polymorphism among 30 samples of H. cordata; 43 (86%) produced fragments of expected size, suggesting that the unigenes were suitable for specific primer design and of high quality, and the SSR marker could be widely used in marker-assisted selection and molecular breeding of H. cordata in the future. CONCLUSIONS: This is the first application of Illumina paired-end sequencing technology to investigate the whole transcriptome of H. cordata and to assemble RNA-seq reads without a reference genome. These data should help researchers investigating the evolution and biological processes of this species. The SSR markers developed can be used for construction of high-resolution genetic linkage maps and for gene-based association analyses in H. cordata. This work will enable future functional genomic research and research into the distinctive active constituents of this genus.

Saururaceae, a member of the paleoherbs, is an ancient family with six species in four genera, Anemopsis, Gymnotheca, Houttuynia and Saururus [1] . Houttuynia cordata Thunb. (Yuxingcao in Chinese) is the only species in the genus Houttuynia [2, 3] . It is distributed mainly in the central, southeastern and southwestern regions of China, and extends to Japan, Korea and Southeast Asia, where it grows in moist, shady places [4] . H. cordata is an important traditional medical herb native to China and other Asian countries [5, 6] . It plays a unique role in improving the immune system of patients with severe acute respiratory syndrome (SARS) [7, 8] . Extracts of H. cordata have diverse pharmacological effects including anticestodal [9] , antibacterial [10, 11] , antiviral [12] [13] [14] [15] , anticancer [16, 17] , antioxidant [18, 19] , antiallergenic [20, 21] , anti-inflammatory [22] [23] [24] , antimutagenic [18] and anti-obesity [25] activities. H. cordata is also consumed as a vegetable in China for its special aroma. Although H. cordata is of high medicinal and nutritional value, there are no genomic resources for this non-model genus. This lack of genomic information has become a limitation for extensive and intensive research on this important traditional medical herb.

Previous studies on this plant have mainly focused on cultivation techniques [26, 27] , its physiological and biochemical properties [28, 29] , its genetic relationships and the diversity among H. cordata germplasm collections from different places [4, 30, 31] , and its pharmacological effects [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] . To date, few gene sequences or novel genes have been reported on this species, although much effort has been devoted to cloning key genes. RNA-Seq, which is based on next generation sequencing, is a high throughput technology that has great advantages in examining the fine structure of a transcriptome [32] . When no genome sequence is available, transcriptome sequencing provides an effective way to obtain large amounts of sequence data [33] . RNA-Seq has been widely used in many organisms to obtain mass sequence data for transcriptional analysis, gene discovery and molecular marker development [34] [35] [36] . The genetic relationships and diversity among H. cordata germplasm collections have been investigated mostly using AFLP [37] , ISSR [4] , PCR-RFLP [30] and RAPD markers [31] . No simple sequence repeat (SSR) markers have been reported in H. cordata. Compared with other types of molecular markers, SSR markers have many advantages, such as simplicity, effectiveness, abundance, hypervariability, reproducibility, codominant inheritance, and extensive genomic coverage [38] . Because of the lack of effective molecular markers, marker-assisted selection and molecular breeding of H. cordata has lagged behind other medicinal plants such as Panax notoginseng (Burkill) F.H.Chen [39] , Gastrodieae elata Blume [40] , and Glycyrrhiza uralensis Fisch. [41] . Thus, a rapid, low-cost and effective approach is required to develop SSRs molecular markers for H. cordata.

In this study, we applied the next-generation massively parallel sequencing technique (Illumina HiSeq 2000) to the sequencing and analysis of the complete H. cordata transcriptome for the first time. We sampled the pooled transcriptomes of flower, leaf, stem and rhizome tissues of H. cordata and used Illumina paired-end sequencing technology to generate a large-scale EST database and develop a set of SSR markers. These results provide a very useful genomic resource for research on H. cordata in the future.

To obtain a global overview of the H. cordata transcriptome and gene activity at the nucleotide resolution, RNA was extracted from four different H. cordata tissues including the rhizome, stem, leaf and flower, and mixed at equivalent concentrations. Sequencing was performed on an Illumina HiSeq2000 genome analyzer. Each sequencing pass can yield two 90-bp independent reads from either end of a DNA fragment. In all, 56,668,324 sequence reads were generated, of which 51,973,070 were of acceptable quality after cleaning the low-quality reads (Table 1) , then we used the Trinity short reads assembly program to assemble the reads for non-redundant consensus [42] 

Several complementary methods were used to annotate the assembled sequences. First, the assembled sequences were queried against the National Center for Biotechnology Information (NCBI) nonredundant protein (Nr) and Swiss-Prot protein databases using BLASTx to search for similar sequences (E-value ,10 25 ). Of the 63,954 assembled sequences, 39,982 (62.52%) showed homology to sequences in the Nr database (Table S1 ), while 26,122 (40.84%) unigenes had homology to proteins in the Swiss-Prot database. In addition, 99.06% of the unigenes over 1,000 bp in length showed homologous matches, whereas only 29.28% of the unigenes shorter than 300 bp showed matches (Fig. 2) . The 2,624 unigenes that had no matches in either the Nr or Swiss-Prot databases were subjected to gene prediction analysis using ESTScan (Ver-sion3.0.2) [43] . In total, 42,785 unigenes were detected by homology analysis using the Nr and Swiss-Prot databases or ESTScan prediction.

The E-value distribution of the top hits in the Nr database showed that 62.2% of the mapped sequences had strong homology (E-value ,10 230 ) and 50.8% had very strong homology (E-value ,10 245 ) to available plant sequences, whereas 37.8% of the homologous sequences had E-values in the range 10 25 to 10 230 (Fig. 3A) . The similarity distribution of sequences database showed that 3,880 (9.70%), 12,513 (31.30%), 16 biaceae), Populus trichocarpa (Salicaceae), Glycine max, Medicago truncatula (Leguminosae), Oryza sativa Japonica Group, and Sorghum bicolor (Gramineae), respectively. And 18.35% of the distinct sequences had matches to sequences from 'other' species. (Fig. 3C ). Based on Nr annotation, 30,131 unigenes were assigned gene ontology (GO) terms. The sequences that belonged to the biological process, cellular component, and molecular function clusters were categorized into 46 functional groups (Fig. 4) . 'Cellular processes' and 'metabolic processes', 'cell' and 'cell part', 'antioxidant activity' and 'binding' were the dominant groups among the three main categories (biological process, cellular component and molecular function), respectively. However, we did not find any genes in the clusters 'carbon utilization', 'locomotion', 'nitrogen utilization', 'sulfur utilization', 'viral reproduction', 'extracellular matrix', 'extracellular matrix part', 'extracellular region part', 'channel regulator activity', 'metallochaperone activity', 'protein tag' or 'translation regulator activity' (Fig. 4) .

In addition, all unigenes were subjected to a search against the Cluster of Orthologous Groups (COG) database for functional prediction and classification. Overall, 15,363 of the 39,982 sequences showing Nr hits were assigned COG classifications. The COG-annotated putative proteins were functionally classified into at least 25 molecular families including cellular structure, biochemistry metabolism, molecular processing, and signal transduction (Fig. 5) .

The cluster for 'general function prediction only' (5,362, 13.58%) represented the largest group, followed by transcription (4,124, 10.44%), replication, recombination and repair (3, 242, 8 .21%), posttranslational modification, protein turnover and chaperones (3,014, 7.63%), function unknown (2,787, 7.06%), signal transduction mechanisms (2,614, 6.62%), translation, ribosomal structure and biogenesis (2,603, 6.59%), cell cycle control, cell division, chromosome partitioning (2,446, 6.19%), carbohydrate transport and metabolism (2,340, 5.93%), cell wall/ membrane/envelope biogenesis (2,263, 5.73%), intracellular trafficking, secretion, and vesicular transport (1,244,3.15%), amino acid transport and metabolism (1,177,2.98%), secondary metabolites biosynthesis, transport and catabolism (1,109, 2.81%), inorganic ion transport and metabolism (909, 2.30%), lipid transport and metabolism (887, 2.25%), energy production and conversion (823, 2.08%), coenzyme transport and metabolism (501, 1.27%), defense mechanisms (447, 1.13%), cytoskeleton cell (398, 1.00%), motility (385, 0.98%), nucleotide transport and metabolism (298, 0.75%), chromatin structure and dynamics (270, 0.68%) and RNA processing and modification (222, 0.56%), whereas only a few unigenes were assigned to nuclear structure and extracellular structure (13 and 11 unigenes, respectively). The COG function classification of H. cordata unigenes is shown in 

A total of 24,434 assembled sequences were associated with 128 predicted KEGG metabolic pathways. The number of sequences ranged from 3 to 6,718. The top 20 pathways with the greatest number of sequences are shown in Table 2 . The greatest number of transcripts was found in the metabolic pathways. The top 10 metabolic pathways were: glycerophospholipid metabolism (1,974), ether lipid metabolism (1,853), starch and sucrose metabolism (869), purine metabolism (697), pyrimidine metabolism (637), phenylpropanoid biosynthesis (381), oxidative phosphorylation (298), amino sugar and nucleotide sugar metabolism (Table S2 ).

For further assessment of the assembly quality and development of new molecular markers, all 63,954 unigenes generated in this study were used to mine potential microsatellites, which were defined as di-to hexa-nucleotide SSRs with a minimum of four repetitions for all motifs. The SSRs that were only located in one single read had been eliminated. Using the MIcroSAtellite (MISA, http://pgrc.ipk-gatersleben.de/misa/) tool, a total of 4,800 potential SSRs were identified in 4,413 unigenes, 357 of which contained more than one SSR; 164 SSRs were present in compound form ( Table 3) .

The frequency, type and distribution of the 4,800 potential SSRs were also analyzed in this study. The compilation of all SSRs revealed that, on average, one SSR could be found every 9.04 kb in the unigenes. Among the 4,800 SSRs, tri-nucleotide repeat motifs were the most abundant type (1,994, 41.54%), followed by mono-(1,313, 27.35%), di-(1,278, 26.63%), hexa-(110, 2.29%), penta-(66, 1.38%) and tetra-nucleotide (39, 0.81%) repeat motifs. The mono-to hexa-nucleotide motifs were further analyzed for SSR repeat numbers (Table 4) .

SSR length was mostly distributed in the 12 to 20 bp range, accounting for 87.04% (4035 SSRs) of the total SSRs, followed by 21-99 bp (562 SSRs, 12.12%). There were 39 SSRs longer than 100 bp.

Within the SSRs, 119 motif sequence types were identified, of which mono-, di-, tri-, tetra-, penta-and hexa-nucleotide repeats had 2, 4, 10, 11, 31 and 58 types, respectively. The most abundant motif detected in the SSRs was the A/T mono-nucleotide repeat (1,296, 27. 

Fifty primer pairs (designated HM_1-HM_50) were randomly selected from the microsatellites, excluding mono-nucleotide repeats motif, to evaluate their applicability and the polymorphism across 30 individuals of H. cordata (Table S3) Table 5 ).

Next generation sequencing technologies provide a low cost, labor saving and rapid means of transcriptome sequencing and characterization [44] , which enables various functional genomic studies on an organism. Although the 454 Life Sciences (Roche) technology is often used for transcriptome analysis of non-model organisms, it is more expensive than the Illumina technology [45] . De novo assembly of short reads without a known reference is considered difficult [46] , but de novo assembly of transcriptomes using short reads has received attention [47] . In this study, we demonstrated a strategy for de novo assembly of a transcriptome using short reads for a non-model medicinal plant, H. cordata, for which sequence data is very limited in the public databases at present. We showed that assembly program parameters and sequence quality have a significant effect on the assembly output.

Although the length of contigs were often less than 500 bp, the Illumina sequencing solution was reliable. such as the average contig size of sesame was less than 200 bp [48] , whitefly was only 40 bp [49] , sweetpotato was 202 bp [50] . Compared with these reports, the assembled contigs in this study was quite long (253 bp). This suggested that the coverage was relatively high. Greater N50 and average lengths are considered indicative of a better assembly. Here, the N50 length of the unigenes was 1,051 bp and the average length was 679 bp, which suggests that the relatively short reads from Illumina paired-end sequencing for this non-model organism have been effectively and accurately assembled. Illumina sequencing yielded 56.67 million paired-end reads for H. cordata. The 63,954 unigenes produced here may be useful for further research into H. cordata functional genomics. Of the H. cordata unigenes, 39,982 (62.52%) showed homology to sequences in the Nr database. Comparatively, in Epimedium sagittatum [51] , whitefly [49] , sweet potato [50] and sesame [48] , only 38.50%, 16.20%, 46.21% and 54.03% of the unigenes, respectively, had homologs in the Nr database. The average unigene length in our database was 679 bp, compared with 246, 266, 581 and 629 bp, respectively, in the four studies mentioned above. The higher percentage of hits found in this study was partially a result of the increased number of long sequences in our unigene database; the results for whitefly [49] and sesame [48] support this conclusion. Homologs in other species were not found for 18.3% of the unique sequences. Specifically, only 29.28% of the unigenes shorter than 300 bp showed matches, meaning that 70.72% produced no hits (Fig. 3) . These shorter sequences may lack a characterized protein domain, or they may contain a known protein domain but the query sequence is too short to show sequence matches, resulting in false-negative results. Additionally, little genomic and transcriptomic information is currently available for H. cordata, and consequently, many H. cordata lineage-specific genes might not be included in current databases.

Both gene annotation and KEGG pathway analyses are useful for predicting potential genes and their functions at a wholetranscriptome level. In the H. cordata transcriptome, the predominant gene clusters are involved in the cellular process and metabolic process categories of the biological process GO domain, the cell and cell part categories of the cellular component domain, and antioxidant activity and binding categories of the molecular function domain. Similar results were found in sesame [51] and whitefly [49] . However, in the chickpea transcriptome, the sequences were found to be mainly involved in protein metabolism (biological process) and transferase activity (molecular function) [52] . This suggests remarkable differences among different species of plants.

KEGG analysis showed that 24,434 sequences were involved in 128 known metabolic or signaling pathways, including endocytosis, plant hormone signal transduction and plant-pathogen interaction. H. cordata is one of the most important medicinal plants and is rich in secondary metabolites, which makes it a very important target for genomic studies. In this study, 2,448 (10.02%) sequences of H. cordata were associated with biosynthesis of secondary metabolites (Table 2) . These results may be useful for further investigation of gene function in the future.

The large number of sequences generated for H. cordata in this study for the first time provides valuable sequence information at the transcriptomic level for screening of novel functional genes, or for investigation of molecular mechanisms.

SSR markers play an important role in genetic diversity research, population genetics, linkage mapping, comparative genomics, and association analysis [38, 53] . Previously, genetic diversity analysis of H. cordata germplasm was restricted to AFLP [37] , ISSR [4] , PCR-RFLP [30] and RAPD markers [31] . One of the main reasons for this was the lack of a genome sequence or transcriptome information for H. cordata. Our results have resolved this problem and enabled development of SSR markers for this species. In the present study, 4,800 perfect microsatellites exceeding 12 bp were identified from the H. cordata dataset, and 119 motif sequence types were identified. If mono-nucleotide repeats were excluded, di-nucleotide repeats were the most abundant type, followed by tri-nucleotide repeats, which is consistent with previous reports [54] [55] [56] [57] . The fact that the most abundant di-and tri-nucleotide motifs were AG/TC and AAG/ TTC, respectively, was also coincident with previous reports on other species of plants [54] [55] [56] .

In this study, 48 (96%) of the primer pairs designed from the unigenes successfully yielded high-quality amplicons. These results suggested that the unigenes were suitable for specific primer design, that the assembled unigenes were of high quality, and that the SSRs identified from our dataset could be useful in the future. Five primer pairs produced products that deviated from the expected size, which might have been caused by the presence of introns [51, 58, 59] , large insertions or repeat number variations, or a lack of specificity [51] . The failure of two primer pairs to produce amplicons might have been because of the primer(s) being located across splice sites, large introns, chimeric primer(s), poorquality sequences, or assembly errors [51, 59] .

The 43 primer pairs in our dataset were used to investigate polymorphisms in 30 individuals of H. cordata from 15 populations located across the natural distribution of the species in 13 provinces of China. The results indicated that the level of polymorphism was relatively high, which was also coincident with previous reports using ISSR [4] , RAPD [31] and RAMP markers [60] . Since we identified other SSRs in our dataset, more PCR primers could be developed that would be very useful in germplasm polymorphism assessment, quantitative trait loci mapping [61] , comparative genomics [62] , functional genomics and proteomics studies.

In this study, we have analyzed the transcriptome of H. cordata using high-throughput Illumina paired-end sequencing. We obtained 39,982 sequences and demonstrated some important features of the H. cordata transcriptome through gene annotation and KEGG pathway analysis. In addition, we identified reliable genetic markers in the form of 4,800 SSRs. Fifty primer pairs were randomly selected to detect polymorphism among 30 H. cordata accessions, and 43 (86%) of these primer pairs successfully amplified fragments, revealing abundant polymorphisms. The SSR markers developed in this study can be used for construction of high-resolution genetic linkage maps and to perform gene-based association analyses in H. cordata. This is the first application of Illumina paired-end sequencing technology to investigate the whole transcriptome of H. cordata and to assemble RNA-seq reads without a reference genome. This study will provide useful resources and markers for functional genomics and proteomics research on H. cordata in future.

H. cordata, whose seeds came from national forest park of Zhong Po Mountain in Huaihua City of Hunan Province, was grown in the experimental station of the Department of Life Sciences, Huaihua University, Huaihua, China. The individuals have 18 chromosomes and diploid karyotype (Fig.7) . Flower, leaf, stem and rhizome tissues were harvested 14 weeks post planting,because H. cordata was planted in spring, and their flowers began to open after 14 weeks. All of the samples were immediately frozen in liquid nitrogen and stored at 280uC until use.

Total RNA was isolated using the TRIzol reagent according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA). Table 5 . Cont.

Forward primer Reverse primer No. of alleles It was treated with RNase-free DNase I (Worthington, Lakewood, CO, USA) for 30 min at 37uC to remove any residual DNA. The quality and quantity of each RNA sample was determined by biophotometer (Eppendorf, Germany). Only those RNA samples with an A 260 :A 280 ratio from 1.9 to 2.1 and an A 260 :A 230 ratio from 2.0 to 2.5 were used for the analysis. A total of 20 mg of RNA was equally pooled from the four different tissues for cDNA library preparation.

Beads with oligo(dT) were used to isolate poly(A) mRNA after total RNA was collected from the samples. Fragmentation buffer was added to disrupt the mRNA into short fragments. Taking these short fragments as templates, random hexamer primers were used to synthesize first-strand cDNA. Second-strand cDNA was then synthesized using buffer, dNTPs, RNase H and DNA polymerase I. The short fragments were purified with a QiaQuick PCR extraction kit and resolved with EB buffer for end reparation and A tailing. The short fragments were then connected with sequencing adapters. After agarose gel electrophoresis, suitable fragments were selected for PCR amplification as templates. The cDNA library was sequenced on an Illumina HiSeq2000 sequencing platform.

The raw reads were cleaned by removing adapter and lowquality sequences, because sequencing errors can create difficulties for the short-read assembly algorithm. We therefore carried out a stringent filtering process. Firstly, we discarded all reads with adaptor contamination. Secondly, we ruled out low-quality reads with ambiguous sequences ''N''. Thirdly, the reads with more than 10% Q,20 bases were also removed. De novo transcriptome assembly was carried out with the short reads assembly program in the Trinity software(Release-20120608) [42] . Contigs were created by combining reads that had a certain length of overlap. The reads were then mapped back to the contigs; with paired-end reads we were able to detect contigs from the same transcript as well as the distances between the contigs. The contigs were connected using the Trinity software to get sequences that could not be extended at either end. Such sequences were defined as unigenes. These unigenes were further processed by sequence splicing and redundancy removal using the TGICL software(Version 2.1) [63] to acquire non-redundant unigenes that were as long as possible. After gene family clustering, the unigenes could be divided to two classes. One was clusters, including unigenes that were .70% similar to each other; the other was singletons. In the final step, the sequence direction of the unigenes was determined.

The unigenes were first aligned to sequences in the NCBI Nr and Swiss-Prot protein databases with an E-value ,10 25 using BLASTx. Unigenes that did not have homologs in the databases were scanned using ESTScan (Version3.0.2) [43] . Blast2GO(Version 2.5.0) [64] was used to obtain GO annotations for the unigenes based on BLASTx hits against the NCBI Nr database with an E-value threshold of ,10 25 . WEGO [65] was used for GO functional classification of all unigenes and to plot the distribution of the H. cordata gene functions. The unigene sequences were also aligned to the COG database to predict and classify their functions. Pathway assignments were carried out based on the KEGG database [66] , which contains a systematic analysis of inner-cell metabolic pathways and the functions of gene products.

A microsatellite program (MISA; http://pgrc.ipkgatersleben. de/misa/) was used to identify and localize microsatellite motifs. We searched for all types of simple sequence repeats (SSRs) from mononucleotide to hexanucleotide using the following parameters: at least 12 repeats for mono-, six repeats for di-, five repeats for triand tetra-, and four repeats for penta-and hexa-nucleotide simple repeats. Primer pairs were designed using the software Primer 3-2.2.2. The major parameters for primer design were set as follows: primer length of 18-25 bases (optimal 21 bases), PCR product size of 80-200 bp (optimal 100-150 bp), GC content of 40-60% (optimal 50%), and DNA melting temperature of 57-64uC (optimal annealing temperature 55-59uC).

Thirty individuals of H. cordata from 15 populations located across the natural distribution of the species in 13 provinces of China (Table S3) were selected for polymorphism investigation with the SSRs. Leaf samples were collected, dried and preserved in silica gel until DNA extraction. Genomic DNA was extracted from the leaves of each individual using the CTAB protocol [67] , dissolved in double distilled water, and quantified using agarose gel electrophoresis. The DNA concentration was calculated according to DNA standards. PCR amplification was performed in 16 mL reaction mixtures. Each reaction contained 0.2 mL Taq DNA polymerase (0.5 U/mL), 2.5 mL PCR buffer, 1.5 mL MgCl 2 (25 mmol/L), 0.5 mL dNTPs (2.5 mmol/L), 0.4 mL each primer (10 pmol/L), 2.0 mL template DNA (50 ng/mL), and 8.5 mL sterilized H 2 O. The temperature profiles were: initial denaturation at 94uC for 3 min, 35 cycles of denaturation at 94uC for 30 s, annealing at the optimal temperature of each primer pair for 30 s, and extension at 72uC for 45 s. Final extension was at 72uC for 5 min, and then samples were held at 4uC. After PCR amplification, 6 mL aliquots of the amplified PCR products were loaded onto an 8% polyacrylamide gel. After 3-4 h of electrophoresis (250 V), the gels were stained with silver nitrate (silver staining). Perfectly amplified loci were tested for polymorphism by genotyping 30 individuals of H. cordata. The genetic diversity and mean allele number were calculated using Popgene version 1.32 [68] . Polymorphism information content (PIC) was derived according to the following formula:

where n is the number of alleles at one locus; Pi and Pj are the frequencies of the ith and jth alleles at one locus; j = i+1.

Table S1

The statistics and annotations of unigenes.

(XLS) 

On the evolution and distribution in Saururaceae

Florae Reipublicae Popularis Sinicae

Genetic Variation and Population Differentiation in a Medical Herb Houttuynia cordata in China Revealed by Inter-Simple Sequence Repeats (ISSRs)

Adverse events to Houttuynia injection: A systematic review

Genetic diversity among the germplasm resources of the genus Houttuynia Thunb. in China based on RAMP markers

Anti-inflammatory effect of Houttuynia cordata injection

Immunomodulatory and anti-SARS activities of Houttuynia cordata

Temjenmongla Anticestodal activity of Houttuynia cordata leaf extract against Hymenolepis diminuta in experimentally infected rats

Variation in chemical composition and antibacterial activities of essential oils from two species of Houttuynia cordata Thunb

Biological and antibacterial activities of the natural herb Houttuynia cordata water extract against the intracellular bacterial pathogen salmonella within the RAW 264.7 macrophage

Virucidal effects of the steam distillate from Houttuynia cordata and its components on HSV-1, influenza virus, and HIV

Houttuynia cordata blocks HSV infection through inhibition of NF-kB activation

vitro and in vivo effects of Houttuynia cordata on infectious bronchitis virus

Pharmacodynamic experiment of the antivirus effect of Houttuynia cordata injection on influenza virus in mice

Houttuyninum, an active constituent of Chinese herbal medicine, inhibits phosphorylation of HER2/neu receptor tyrosine kinase and the tumor growth of HER2/neuoverexpressing cancer cells

Houttuynia cordata Thunb. extract inhibits cell growth and induces apoptosis in human primary colorectal cancer cells

A study of the antioxidative and antimutagenic effects of Houttuynia cordata Thunb. using an oxidized frying oil-fed model

Antioxidative effects of polyphenols in leaves of Houttuynia cordata on protein fragmentation by copper-hydrogen peroxide in vitro

Inhibitory effects of Houttuynia cordata water extracts on anaphylactic reaction and mast cell activation

Suppressive effects of Houttuynia cordata Thunb. (Saururaceae) extract on Th2 immune response

Houttuynia cordata, a novel and selective COX-2 inhibitor with anti-inflammatory activity

Anti-inflammatory effects of Houttuynia cordata supercritical extract in carrageenan-air pouch inflammation model

Anti-inflammatory effects of a Houttuynia cordata supercritical extract

Water extract of Houttuynia cordata Thunb. leaves exerts anti-obesity effects by inhibiting fatty acid and glycerol absorption

Study on Highly Effective Cultivation Techniques of Houttuynia cordata Thunb in Cold Alpine Region

Study on germplasm and planting standard operation procedures of Houttuynia cordata Thunb

Moisture gradient affects the growth and physiological characteristics of Houttuynia cordata

Phytohistochemistry on Developing Anthers of Houttuynia cordata

PCR-RFLP analysis of cpDNA and mtDNA in the genus Houttuynia in some areas of China

RAPD analysis on the germplasm resources of herba Houttuynia

Transcriptomic analysis of Chinese bayberry (Myrica rubra) fruit development and ripening using RNA-Seq

Transcriptome Analysis of Silver Carp(Hypophthalmichthys molitrix) by Paired-End RNA Sequencing

Critical assessment of assembly strategies for non-model species mRNA-Seq data and application of next-generation sequencing to the comparison of C(3) and C(4) species

Gene discovery using massively parallel pyrosequencing to develop ESTs for the flesh fly Sarcophaga crassipalpis

Deep sequencingbased transcriptome profiling analysis of bacteria-challenged Lateolabrax japonicus reveals insight into the immune-relevant genes in marine fish

Analysis on genetic diversity of different geographical populations of Houttuynia cordata

The comparison of RFLP, RAPD, AFLP and SSR (microsatellite) markers for germplasm analysis

Breeding strains of Panax notoginseng by using EST-SSR markers

Analysis in different gemplasms of Gastrodieae elata Bl.based on AFLP and SSR markers

Analysis of SSR Information in EST Resource of Glycyrrhiza Uralensis

Fulllength transcriptome assembly from RNA-Seq data without a reference genome

ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences

Applications of new sequencing technologies for transcriptome analysis

Rapid transcriptome characterization for a non-model organism using 454 pyrosequencing

Next-generation sequencing transforms today's biology

Benchmarking next-generation transcriptome sequencing for functional and evolutionary genomics

Characterization of the sesame (Sesamum indicum L.) global transcriptome using Illumina pairedend sequencing and development of EST-SSR markers

Deseea rncho arvtiocle characterization of a whitefly transcriptome and analysis of its gene expression during development

De novo assembly and characterization of root transcriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato (Ipomoea batatas)

Development of a EST dataset and characterization of EST-SSRs in a traditional Chinese medicinal plant, Epimedium sagittatum

De Novo Assembly of Chickpea Transcriptome Using Short Reads for Gene Discovery and Marker Identification

The power and promise of population genomics: from genotyping to genome typing

De novo assembly and characterization of bark transcriptome using Illumina sequencing and development of EST-SSR markers in rubber tree

Transcriptome sequencing of Hevea brasiliensis for development of microsatellite markers and construction of a genetic linkage map

Development, characterization and cross-species/genera transferability of EST-SSR markers for rubber tree (Hevea brasiliensis)

Development and application of EST-SSR markers in Hevea brasiliensis Muell

Tall fescue EST-SSR markers with transferability across several grass species

Genic microsatellite markers in plants: features and applications

Genetic diversity among the germplasm resources of the genus Houttuynia Thunb. in China based on RAMP markers

QTL mapping of ten agronomic traits on the soybean (Glycine max L. Merr.) genetic map and their association with EST markers

The genome sequence of Caenorhabditis briggsae: a platform for comparative genomics

TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets

Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research

WEGO: a web tool for plotting GO annotations

KEGG: kyoto encyclopedia of genes and genomes

DNA Protocols for Plants CTAB Total DNA Isolation

Population genetic analysis of co-dominant and dominant markers and quantitative traits

We thank the anonymous referees and the editor for their comments and suggestions that helped improve the manuscript.