{'fraction_non_alphanumeric': 0.032639151836425594, 'fraction_numerical': 0.03949261643316925, 'mean_word_length': 4.7253414264036415, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transcriptome Sequencing of Hevea brasiliensis for Development of Microsatellite Markers and Construction of a Genetic Linkage Map

To obtain more information on the Hevea brasiliensis genome, we sequenced the transcriptome from the vegetative shoot apex yielding 2 311 497 reads. Clustering and assembly of the reads produced a total of 113 313 unique sequences, comprising 28 387 isotigs and 84 926 singletons. Also, 17 819 expressed sequence tag (EST)-simple sequence repeats (SSRs) were identified from the data set. To demonstrate the use of this EST resource for marker development, primers were designed for 430 of the EST-SSRs. Three hundred and twenty-three primer pairs were amplifiable in H. brasiliensis clones. Polymorphic information content values of selected 47 SSRs among 20 H. brasiliensis clones ranged from 0.13 to 0.71, with an average of 0.51. A dendrogram of genetic similarities between the 20 H. brasiliensis clones using these 47 EST-SSRs suggested two distinct groups that correlated well with clone pedigree. These novel EST-SSRs together with the published SSRs were used for the construction of an integrated parental linkage map of H. brasiliensis based on 81 lines of an F1 mapping population. The map consisted of 97 loci, consisting of 37 novel EST-SSRs and 60 published SSRs, distributed on 23 linkage groups and covered 842.9 cM with a mean interval of 11.9 cM and ∼4 loci per linkage group. Although the numbers of linkage groups exceed the haploid number (18), but with several common markers between homologous linkage groups with the previous map indicated that the F1 map in this study is appropriate for further study in marker-assisted selection.


Introduction
Hevea brasiliensis, commonly known as rubber tree, is almost the sole source of natural rubber production. Natural rubber has a wide range of industrial applications and is under increasing global demand. Hevea brasiliensis is a perennial cross-pollinating and monoecious plant that belongs to the Euphorbiaceae family. The observation of tetravalents during meiosis has lead to the conclusion that H. brasiliensis is a stabilized amphidiploid (2n ¼ 4x ¼ 36). 1 However, the pattern of marker ratios segregating in a population of over 100 trees suggests that H. brasiliensis behaves as a diploid (2n ¼ 36). 2 Several research groups have developed molecular markers to study the genetic diversity of H. brasiliensis, 3 -8 including isozymes, 9 restriction fragment length polymorphisms (RFLPs), 10 amplified fragment length polymorphisms (AFLPs), 2 microsatellites (or simple sequence repeats, SSRs) 2,8,11,12 and expressed sequence tag (EST)-SSRs. 3 These markers have also been used to construct the linkage maps 2 and quantitative trait loci (QTL) maps. 13 -16 Lespinasse et al. 14 produced the first rubber tree linkage map containing 301 RFLPs, 388 AFLPs, 18 SSRs and 10 isozymes, which was later used to identify the QTL variants conferring resistance to the South American leaf blight. 13 -15 Recently, Le Guen et al. 16 constructed the linkage maps based on SSR and AFLP markers and were able to identify the QTL conferring resistance to Microcyclus ulei.
In many organisms, ESTs have been useful for the annotation of genes during genome sequencing efforts, 17 for comparative genome studies 18 and for the production of a genetic linkage map. 19 To date, there are only 12 365 ESTs from H. brasiliensis in GenBank, restricting the quality of research that can be performed on this important plant species. Previous transcriptome studies of H. brasiliensis have been limited in range, focusing mainly on latex in order to gain insight into the rubber biosynthesis pathways. 20 -22 Of the available H. brasiliensis ESTs, 11 256 ESTs are from latex, 1091 ESTs are from bark and 18 ESTs are from leaves. In addition to gene discovery, EST resources enable the identification of markers such as EST-SSRs and single-nucleotide polymorphism. Since these markers are directly linked to functional genes, they are useful for assessing genetic diversity and mapping phenotypic traits. Feng et al. 3 identified 799 SSRs in 10 829 ESTs available in the GenBank database and carried out the genetic diversity assessment of H. brasiliensis using 87 EST-SSR markers. The result provided evidence for cross-taxa transferability and indicated moderate polymorphisms of EST-SSR markers in Hevea species. 3 However, additional markers are desirable to enable quality research into the genetic basis of commercially relevant traits that can be used in marker-assisted breeding programs.
Genomic and transcriptomic resources for H. brasiliensis can greatly benefit from the application of the recent high-throughput sequencing technology, such as the 454 pyrosequencer, 23 which has been instrumental in the development of genetic databases for several economical crops. 24 -27 The purpose of the present study is therefore to sequence the transcriptome of the shoot apical tissue, which is a highly dynamic structure, to discover genes, expand the EST database and develop EST-SSR markers that can be used for assessing genetic diversity, constructing linkage maps and identifying traits of commercial interest.

Plant materials
The shoot apical meristem (SAM; 1-cm long from the vegetative shoot apex) of H. brasiliensis (clone RRIM600) was collected for RNA extraction from an experimental field at the Rubber Research Institute of Thailand, Ministry of Agriculture and Cooperatives, Thailand. The sample was immediately frozen in liquid nitrogen and stored at 2808C until RNA extraction. For the analysis of SSR markers, leaf samples from 20 clones of H. brasiliensis, 2 accessions of Manihot esculenta, 3 accessions of Jatropha curcas and 1 accession of Jatropha gossypifolia (Supplementary data, File S1) were collected and DNA was extracted using a DNeasy Plant Mini Kit (Qiagen). For genetic linkage map construction, 81 samples of an F 1 mapping population were developed from a cross between RRIM600 as a female parent and RRII105 as a male parent. The plants were grown at Chachoengsao Rubber Research Center Office of Agriculture Research and Development, Department of Agriculture, Ministry of Agriculture and Cooperatives, Thailand. DNA samples were extracted using a DNeasy Plant Mini Kit (Qiagen). The concentration of each sample was calculated from the OD measurement using a nanodrop ND1000 (NanoDrop Technologies).

cDNA library preparation and sequencing
Total RNA was extracted using Concert TM Plant RNA Reagent (Invitrogen). Two hundred nanograms of the poly-A mRNA sample was isolated using an Absolutely mRNA Purification Kit (Stratagene) and fragmented in 10Â fragmentation buffer (0.1 M ZnCl 2 , 0.1 M Tris -HCl, pH 7.0) at 708C for 30 s. The reaction was stopped by adding 2 ml of 0.5 M EDTA and 28 ml of 10 mM Tris -HCl, pH 7.5. The mRNA sample was cleaned using Agencourt RNAClean reagent (Beckman Coulter), washed with 200 ml of 70% EtOH, air dried and eluted in 20 ml of 10 mM Tris -HCl, pH 7.5. Fragmented mRNA samples were converted to double-stranded cDNA with the cDNA Synthesis System Kit (Roche Applied Sciences) using random primers and AMV Reverse Transcriptase. A cDNA library for 454 pyrosequencing was prepared according to the October 2008 version of the cDNA Rapid Library Preparation protocol (Roche Applied Sciences). The cDNA library was amplified in emulsion PCR and subject to pyrosequencing on two full picotiter plates of the Genome Sequencer (GS) FLX

Sequence analysis
Poly-A/T(18) and 454-adapter sequences were trimmed off. Sequence reads with low quality (average quality scores ,20), short reads (,100 bp), rRNAs and tRNAs were removed. Sequence assembly was performed using the cDNA option of Newbler 2.5, the de novo sequence assembly software. The cDNA option assembles reads into contigs much like a genomic assembly; however, this option allows/expects contigs to have multiple joining contigs representing alternately spliced genes. These contigs (which essentially represent exons) are then assembled into isotigs (representing processed mRNA) and isotigs utilizing the overlapping subsets of contigs are grouped into isogroups (representing genes, isoforms or gene families). When the maximum number of contigs in an isogroup is exceeded, the contigs are output as un-traversed contigs in the isotigs file; all further mention of isotigs will include these un-traversed contigs. Unique sequences were searched for sequence homology against the Uniprot plant protein database (www.uniprot.org); and reference protein sequences of Manihot, Ricinus, Arabidopsis and Oryza (www. phytozome.net) using the BLASTx program with a cutoff at E-6. 28 The assignment of functionality via gene ontology (GO) was performed using Blast2GO. 29 The MISA-MIcro SAtellite identification tool (http://pgrc.ipk-gatersleben.de/misa/misa.html) was used to search for SSR from the EST data set. For the searches and comparison of microsatellites, SSRs were defined as being mononucleotide repeats (MNRs) !10 repeats and di-(DNRs), tri-(TNRs), tetra-(TTNRs), penta-and hexanucleotide repeats !6 repeats; criteria for composite SSRs was an interval of bases 100. For the purpose of marker evaluation, we increased stringency to reduce the number of candidates. We designed primer pairs overlapping DNRs and TNRs !8, TTNRs !7, pentanucleotide repeats (and more) !6 or containing complex SSRs !30 nucleotides. In some cases, candidate SSRs that passed the criteria suggested by Feng et al. 3 were prioritized based on the presence of motif size polymorphisms in the sequence alignment results.

SSR markers from the previous reports
Primer pairs designed for the amplification of genomic SSR markers from the NCBI database (AY486558.1 -AY486910.1) and previously described 8 were used to construct a linkage map together with novel EST-SSR markers.

SSR analysis
DNA samples were extracted from young leaf tissue using a DNeasy Plant Mini Kit (Qiagen). Primer pairs were designed to amplify SSR regions using PRIMER3. 30 PCR was carried out in a total volume of 10 ml containing 2 ng of DNA template, 1Â Taq buffer, 2 mM MgCl 2 , 0.2 mM dNTPs, 1 U Taq-DNA polymerase (Fermentas) and 0.5 mM each of forward and reverse primers. Amplification was performed in a GeneAmp PCR 9700 thermocycler (Applied Biosystems) programmed as follows: 948C for 2 min followed by 35 cycles of 948C for 30 s, 528C for 30 s, 728C for 1 min and a final extension step at 728C for 10 min. Amplified products were separated on 5% denaturing polyacrylamide gels and visualized by silver staining.

Analysis of polymorphic loci
Twenty accessions of H. brasiliensis, as listed in Supplementary data, File S1, were used for the polymorphism analysis of SSR markers. Details of primer pairs for amplifiable EST-SSR markers are listed in Supplementary data, File S2. Scored data from polymorphic loci were used to calculate the polymorphism information content (PIC). 31 Observed heterozygosity and expected heterozygosity were calculated using the PowerMarker 3.25 software. 32 The cross-taxa transferability of H. brasiliensis SSR loci were evaluated using six other taxa of Euphorbiaceae plants, including two accessions of M. esculenta, three accessions of J. curcas and one J. gossypifolia (Supplementary File S1). The percentage of transferability was calculated for each taxon by dividing the number of successfully amplified SSR loci by the total number of loci analysed. A genetic similarity matrix was prepared for the 20 H. brasiliensis genotypes at 47 EST-SSR loci (Supplementary File S3) using the NTSYSpc 2.2 software. 33 UPGMA (un-weighted pair group method with arithmetic mean) cluster analysis was conducted usizSpc 2.2 software. 33

Linkage map construction
Eighty-one of H. brasiliensis progenies derived from across between RRIM600 and RRII105 were used as mapping population. Genomic DNA of individual samples was used to genotype with informative primers and genotypic data were scored as codominant markers under the cross-pollination model, e.g. ,abxcd., ,efxeg., ,lmxll., ,nnxnp. and ,hkxhk.) with up to four distinguishable alleles as described by Van Ooijen and Voorrips. 34 The integrated parental genetic linkage map was constructed using the double pseudo-testcross strategy by JoinMap 3.0. 34 The Mendelian segregation ration of all markers was evaluated using the chi-square test K. Triwitayakorn et al. 473 (x 2 ) and distorted markers (P 0.1) were excluded. The map was constructed with an LOD score threshold of 3.0 and the mapping parameters were set with a recombination threshold of 0.4, a jump threshold of 5.0 and a minimum LOD score threshold of 1.0. The map distance between two markers in centiMorgan (cM) was calculated using Kosambi's mapping function. 35 The linkage map was drawn using MapChart 2.2. 36

Results and discussion
3.1. Transcriptome sequencing of rubber tree A total of 2 311 497 filtered sequence reads were generated from the vegetative shoot apical tissue with an average read length of 294 bases totalling 676.5 Mb. The majority of the reads were in the range 200 -400 bp (Fig. 1). All reads were deposited in DDBJ Read Archive (ID ¼ DRA000170). There were 191 369 (8.27%) reads with homology to plant tRNAs and rRNAs and 139 838 (6.04%) reads shorter than 100 bp which were removed before sequence assembly. Raw sequencing reads were assembled by Newbler 2.5, the de novo sequence assembly software, 23 currently the most robust software for 454 transcriptome assembly. 37 A total of 28 387 isotigs from 19 152 isogroups and 84 926 singletons were obtained from the assembly. An isogroup theoretically represents a single gene; however, genes with high sequence similarity may be grouped together and therefore isogroups may also represent gene isoforms or gene families. The isotigs had an average length of 1326 bp and the majority of isotigs were between 500 and 1000 bp ( Table 1). The largest isotig (isotig08423) was 9041 bp, which showed sequence similarity to the M. esculenta chloroplast polycistronic transcript psaA-psaB (YP_001718437.1; E-value ¼ 0). The most highly represented isotig (isotig00002) was assembled from 22 641 reads and the highest blast match was the Hevea hydroxynitrile lyase Chain A (AAC49184.1; E-value ¼ 1E-59), reflecting the high level of cyanogenesis activity in the tissue sample. It is well recognized that all living tissues of H. brasiliensis, including seeds, are strongly cyanogenic accumulating high quantities of cyanogenic precursors such as linamarin and lotaustralin. 38 The average GC content of the H. brasiliensis transcriptome generated in this study was 42.16%, which is similar to the GC content of H. brasiliensis sequences in the GenBank EST database (42.18%). The GC content of H. brasiliensis coding sequences is slightly lower than the average GC content of Arabidopsis coding sequences (44.5%) and rice coding sequences (51.5%); 39 but much higher than Arabidopsis intergenic regions (32.9%; http//gi.kuicr.kyoto-u.ac.jp). All unique transcripts were annotated and characterized according to GO using BLAST2GO 29 and the result is available at http://www4a.biotec.or.th/rubber.  In total, 61 625 isotigs and singletons were assigned one or more GO terms. Under the biological process domain, 71 071 assignments were made, with a large proportion of assignments falling into the categories metabolic process (31.07%) and cellular process (29.94%). A total of 60 927 assignments were made to the molecular function domain, with the majority falling into the categories binding activity (46.52%) and catalytic activity (39.82%). This distribution of GO terms is similar to the previous study in the pea SAM transcriptome sequencing. 40 The large number of annotated sequences shows that ESTs generated by high-throughput sequencing are more likely to represent mRNA than ESTs generated by lower-throughput methods, such as previous studies performed on the latex transcriptome. 20 -22 The reason for this is that high-throughput sequencing generates sequence data supported by multiple reads, often representing complete mRNAs. The majority of EST sequences currently in GenBank from Chow et al. 20 are mostly un-annotated singletons with only 904 sequences (26%) that have GO terms assigned. Also, the diversity of the latex transcriptome is limited as Han et al. 22 pointed out that the genes expressed in latex are mainly associated with rubber biosynthesis pathways, defence mechanisms and allergenic proteins. A larger diversity of genes was considered more likely to be found in the vegetative shoot apex tissue than in latex, and that is what has been found here.
A unique characteristic of the shoot apical tissue is the maintenance of the SAM via intercellular communication involving a complex signalling network such as epigenetic control, 41,42 transcriptional gene regulation 43 -45 and hormonal regulation. 46 Key regulatory genes controlling SAM maintenance, such as WUSCHEL (WUS) and SHOOT MERISTEMLESS (STM) genes, were identified as putative full-length cDNAs in this data set. Recent studies revealed that WUS plays the key role in SAM maintenance through the regulatory loop of WUS -CLAVATA (CLV) feedback 41,47,48 and interacts with STM via phytohormone signalling pathways. 49 Furthermore, the regulation of WUS expression is also controlled by auxin signalling, chromatin remodelling and positive and negative transcriptional regulators. 50 Transcripts for a positive transcriptional regulator of WUS, such as APETELA2, 45 SPLAYED and BARD1, which bind to the WUS promoter sequence 43 were found in this data set. Whereas, WUS negative transcriptional regulators which are required for development of floral organs, such as ULTRAPETALA1 51 and HANABA TARANU 52 transcripts, were not detected.
Also identified were genes from the KNOX (Knottedlike homeobox) family such as STM and KNOTTED-LIKE FROM ARABIDOPSIS THALIANA (KNAT) genes which are essential in maintaining the balance between organ primordia growth and stem cell maintenance in the SAM. 53 KNOX transcription factors have roles in suppressing gibberellins (GAs) in the SAM by inhibiting GA-20 oxidase, which is required for GA biosynthesis, and promoting GA-2 oxidase, which inactivates the active GA. 54 Thus, KNOX proteins prevent the accumulation of GA in the central zone of the SAM, consequently preventing the differentiation of stem cells. KNOX proteins also promote cytokinin activity in the SAM central zone stimulating division and maintenance of undifferentiated stem cells. 46 Flanking the SAM, high levels of auxin and GA activities have roles in development and growth of lateral organ primordia. 46 The tissue sample in this study contained both undifferentiated meristem and differentiated organ primordia; therefore, it was likely to identify transcripts of genes involving in many phytohormone biosynthesis and signalling pathways.

Similarity of rubber tree ESTs to other plant
proteins To investigate the efficiency of gene discovery in the H. brasiliensis transcriptome, the isotigs and singletons were searched for homology using BLAST 28    should also be noted that sequence homology analysis by BLAST can be biased by the number and quality of query database and the reference databases.

EST-SSR: distribution and frequencies
A total of 17 819 SSRs were identified in the isotigs and singletons from the H. brasiliensis transcriptome ( Table 2). This represents an average frequency of one EST-SSR in every 3383 bp, which is a lower frequency than the previous report (1 SSR per 2.25 kb) by Feng et al. 3 Among plant species, SSR frequencies range from 1 per 1.5 kb in coffee 57 to 1 per 67 kb in mungbean. 25 The distribution and frequencies of EST-SSRs significantly vary between different studies due to SSR search criteria, the size of the EST database and software tools, 58 so the MISA software was used with the same criteria as the previous study on H. brasiliensis EST-SSRs 3 to allow for direct comparability. From 17 819 SSRs identified, there were 12 682 MNRs, 3458 DNRs, 1593 TNRs, 51 TTNRs and 35 SSRs with pentanucleotide repeats or more. It should be noted that the number of MNRs may not be accurate due to the limitations of the 454 technology in reading long homopolymer sequences. MNRs were described here for the purpose of comparison with previous studies but they were not used for polymorphism analysis. The most common type of DNR was AG/CT which accounted for 64.43% of the repeats, followed by AT/TA (26.75%), AC/GT (8.48%) and GC/CG (0.3%). The most common type of TNR was AAG/CTT (32.53%), followed by AAT/ATT (19.17%) and ACC/ GGT (18.33%). The least frequent DNR and TNR motif types were GC-rich motifs, which were found at only 33 loci. SSRs with GC-rich motif repeats are rare in many plants, such as H. brasiliensis, 3 rice, corn, soybean, 59 wheat, 60 Arabidopsis, apricot, peach 61 and coffee. 57

Polymorphism test in EST markers
Four hundred and thirty EST-SSRs were selected from EST-SSR present in the data set to give a range of repeat units and motif types, and primers were designed to amplify them. Based on the sequence homology search, 16 of 430 (3.7%) SSR primer pairs were mapped to the EST reads with GenBank accession numbers reported for the development of EST-SSR markers by Feng et al. 3 The primary PCR screening in three clones of H. brasiliensis (RRIM600, RRIC110 and BPM24) showed that 323 primer pairs (75.11%) were amplifiable ( primer sequences were listed in Supplementary File S2). We selected 47 primer pairs that flanked long SSR motifs (criteria mentioned above), gave clear PCR bands in the primary amplification screening and represented in all classes of nucleotide repeats for polymorphism evaluation among 20 different H. brasiliensis clones ( Table 3). The number of alleles observed at each locus ranged from two to six, with an average of 3.85. Although H. brasiliensis is believed to be a stabilized amphidiploid, only five SSRs (EHB61, EHB100, EHB109, EHB115 and EHB144) gave more than two alleles supporting the report by Lespinasse et al. 14 that H. brasiliensis behaves as a diploid. The value of expected heterozygosity varied from 0.1349 to 0.7494 with an average of 0.5594, while direct count heterozygosity ranged from 0.1429 to 0.8095 with an average of 0.5076. PIC values ranged from 0.13 to 0.71 with an average of 0.50, which is higher than the previous report (PIC ¼ 0.38) by Feng et al. 3 The higher PIC value in this study was probably due to a larger sample size of H. brasiliensis clones tested (n ¼ 20) than the previous study (n ¼ 12). The average number of alleles and PIC values of EST-SSR markers are lower than those of genomic SSR markers, 8 as expected for functional sequences.
The 47 polymorphic EST-SSR markers were used to evaluate the genetic relatedness among 20 different H. brasiliensis clones which were classified into two groups at the level of genetic similarity 0.44 (Fig. 3). This generally corresponded well to the clone pedigree (Supplementary File S1). The first group (Group I) contains clones: PB260, PB310, RRIM605, PB217, PB5/51 and PB235. The majority of these clones have PB5/51 or PB49 clones as one of the parental lines. Group II contains a mixture of primary clones and cultivated clones from various rubber  2999  1915  1435  1127  813  688  4393  12 682   DNR  991  535  414  307  237  200  123  132  110  94  315  3458   TNR  700  292  192  114  103  45  42  31  14  15  45  1593   TTNR  33  10  5       Since these markers are able to reproduce the relationship that was already known from pedigree information, they can be used to provide reliable genotypic information for clonal identification and the selection of parents in breeding programs. The minimal set of five highly informative microsatellites (EHB85, EHB109, EHB169, EHB177 and EHB178) was able to distinguish each of the H. brasiliensis clones included in this study.
Genotyping of the 93 EST-SSR makers that were amplifiable in H. brasiliensis was performed across genera with two accessions of M. esculenta, three accessions of J. curcas and one accession of J. gossypifolia, all belonging to the same family as H. brasiliensis, Euphorbiaceae. The results showed that 47 of 93 primer pairs (40%) gave successful amplifications in Manihot species. Fourteen primer pairs (15%) and nine primer pairs (10%) were successfully amplified in J. curcas and J. gossypifolia, respectively. The higher rate of cross-transferability in Manihot species compared with that in Jatropha species suggests a closer relationship of Hevea with Manihot than with Jatropha. Six primer pairs, EHB61, EHB63, EHB85, EHB115, EHB116 and EHB156, were able to amplify unique products from all plant taxa tested.

Genetic linkage map
From 323 novel EST-SSR primer pairs amplifiable in H. brasiliensis, 59 primer pairs were polymorphic between RRIM600 and RRII105 parental clones. All of these polymorphic primers were used to genotype with 81 individual F 1 samples together with 98 published SSR markers. Genotypic data were scored and subjected to linkage analysis. Of 157 polymorphic markers, 124 markers (78.9%) revealed the expected Mendelian segregation ratio and used to construct the genetic linkage map. The F 1 map consisted of 97 loci distributed on 23 linkage groups. Of these, 37 loci were novel EST-SSRs. The total map distance covered 842.9 cM with a mean interval of 11.9 cM and the average loci per linkage group were approximately four loci (Fig. 4). The number of linkage groups exceed the expected haploid number of linkage groups (18 linkage groups), suggesting that more markers are required to fill the gap between adjacent markers. Moreover, the comparison between the F 1 map in this study and the map constructed by Le Guen et al. 16 using different parents revealed 25 common markers on nine homologous linkage groups between both maps (Table 4). A total of seven marker intervals of 20 markers showed co-linearity between homologous linkage groups. Some linkage groups in this study (LG7-LG8, LG11-LG15 and LG12-LG13) could be joined together based on the linked markers in Le Guen et al. 16 The common and collinear markers indicated the reliability between different maps. 62 Therefore, the F 1 map of this study is appropriate for further studies in marker-assisted selection.