Genomic sequencing and microsatellite marker development for Boswellia papyrifera, an economically important but threatened tree native to dry tropical forests

Plant material

Boswellia papyrifera is one of the six Boswellia species that grow in various parts of Ethiopia. The B. papyrifera genotype used for Illumina paired-end sequencing was collected from a natural population at Kafa Humera Wuhdet (14.05265N latitude; 37.13078E longitude) in north-west Ethiopia. Young leaves were collected from growing shoot tips of the plant and preserved in silica gel while in the field and during transportation to the laboratory for DNA extraction. A genomic DNA library for Illumina paired-end sequencing was prepared from 4 µg of DNA following the PCR-based gel-free illumina TruSeq DNA sample prep protocol and sequenced as 2 × 100 nt paired-end reads on an Illumina HiSeq at Greenomics, Wageningen UR, Wageningen, the Netherlands.

Plant material for SSR marker development

For testing of the SSR loci a set of 12 genotypes were used. Ten of the genotypes represented populations of B. papyrifera collected from 10 different regions of Ethiopia. Two genotypes of Boswellia pirrotae and B. popovina were included for testing the cross-transferability of the markers to closely related species. The B. pirrotae sample was from the northwestern part of Ethiopia. Boswellia popoviana is endemic to Socotra Island, Yemen, and the dried leaf sample was obtained through the Edinburgh Royal Botanical Garden, UK.

DNA extraction

Total DNA was isolated from silica-dried young leaves following the cetyltrimethylammonium bromide protocol of Fulton et al. (1995). As large amounts of phenolic compounds were expected because of the resin content in the leaves, the protocol was modified by the addition of 2 % pvp-40 to the extraction buffer and 1 % mercaptoethanol to the microprep buffer of Fulton et al. (1995), added immediately before use. The extraction was followed by purification steps using DNeasy (Qiagen, Venlo, The Netherlands) according to Smulders et al. (2010). DNA yield and quality were visually assessed on a 1 % agarose gel.

Sequence filtering

The raw reads were error-corrected using musket (Liu et al. 2012b). This error-corrected set was used for the repeat assembly. Prinseq-lite 0.20.04 (Schmieder and Edwards 2011) was used for quality control and filtering of reads (minimum read length of 50 nt, minimum average base quality of 25, maximum ambiguous nt (N) of 1) after which the data were used for SSR mining. After low complexity trimming (minimum DUST score of 7 for removal of low complexity reads and removal of duplicate reads, also with Prinseq-lite), paired-end reads with overlapping sequences were connected using connecting overlapped pair-end (Liu et al. 2012a) in the full mode. Reads were filtered for chloroplast sequences by mapping the reads against the closest chloroplast genome available, which is one of Citrus sinensis, using bowtie2 (Langmead and Salzberg 2012, settings -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -a).

Repeat analysis

Reads from the highly repeated fraction of the genome were extracted and assembled using RepARK (REPetitive motif detection by Assembly of Repetitive k-mers; Koch et al. 2014). The motifs present in the repetitive contigs were counted and analysed by blastn (e-value 1e-5) against Repbase v19.08 (database of repetitive DNA elements, Jurka et al. 2005).

Assembly and annotation

A de novo draft assembly was created from the filtered reads using SOAPdenovo 2.21 (Li et al. 2010, settings -K 41 -M 3 -d 4). The gaps emerging during the scaffolding process by SOAPdenovo were closed using GapCloser (vs. 1.12). The contigs >1000 bp of the draft assembly were analysed and functionally annotated using Blast2GO (Conesa et al. 2005).

SSR mining and design of primers

Five million of the filtered but not assembled reads were analysed with PAL_FINDER 0.02.03 (Castoe et al. 2012) to identify SSRs using slightly adjusted criteria: at least six contiguous repeat units for dinucleotide repeats, four for tri- and tetranucleotide repeats and three for penta- and hexanucleotide repeats (Castoe et al. (2012) used six units for trunicleotide repeats). Following Castoe et al. (2012) the reads with multiple SSR loci were considered a ‘compound’ repeat if the SSRs had a different repeat motif, but a ‘broken’ repeat if the SSRs had the same motif. Reverse-complement repeat motifs (e.g. TG and CA) and translated or shifted motifs (e.g. TGG, GTG and GGT) were grouped together, so that there were a total of four unique dinucleotide repeats, 10 unique trinucleotide repeats and so on.

A subset of over 70 000 trinucleotide to hexanucleotide repeat-containing reads was used to further screen potentially amplifiable SSR loci (PALs): loci for which PCR primers could be designed. Primer designing followed the default parameters specified in Primer3 (Rozen and Skaletsky 2000). The reads were then screened for differences in lengths of those sequences that contained these primers (as in Vukosavljev et al. 2015). At these loci the sequenced plant may be heterozygous, thus indicating that the locus is polymorphic. These formed the group of potentially polymorphic loci.

SSR loci amplification and analyses of polymorphism

PCRs were performed in a total volume of 10 µL reaction mixture containing 4 µL 2 ng µL⁻¹ DNA, 5 µL MP mix from Qiagen kit, 0.8 µL (2 µM) universal fluorescent-labelled primer and 0.2 µL mix of the forward and reverse primers. The fluorescent labelling method described in Schuelke (2000) was adapted to label the primers for analyses of the PCR products with a laser detection system. For this the forward primers were labelled with a universal M13 sequence (AACAGGTATGACCATGA) at the 5′ end while the reverse primers were tailed with GTTT at their 5′ end according to Brownstein et al. (1996) to reduce stutter bands (both tailing sequences are not shown in the sequences in Table 1). A thermal cycling profile was set at 15 min of initial denaturation at 95 °C, followed by 30 cycles of 30 s denaturation at 94 °C, 45 s annealing at 56 °C and 45 s extension at 72 °C. This was followed by additional eight cycles with 53 °C annealing temperature to facilitate the annealing of the fluorescent dye-labelled M13 primer, and a final extension step of 10 min at 72 °C. After amplification 10 µL water was added. Fluorescently labelled amplicons were resolved on a 4200 or 4300 Licor DNA analyser.

Next-generation sequencing

Genomic DNA of one B. papyrifera individual was sequenced in order to obtain a library to mine for microsatellite repeats. One lane on an Illumina HiSeq produced 143 458 368 raw reads. Based on k-mer counts, the estimated genome size of B. papyrifera was 705 Mb, sequenced at 36× coverage. After error correction and filtering reads for short sequences, sequences with ambiguities (Ns) and low complexity, and excluding redundant sequences, 120 479 203(84 %) paired-end reads and 10 851 777 single-end reads remained.

SSR identification

A search of SSRs in a subset of five million Illumina paired-end reads identified 170 832 reads (3.4 %) containing SSRs. In these reads, a total of 175 607 repeat loci (dinucleotide through hexanucleotide repeats) were identified, which corresponds to one SSR locus per 5.7 kb. Figure 1 shows the frequency of the top-20 repeat motifs. These include all dinucleotide motif repeats (of at least six repeat units long), of which AC and AT repeats were the most abundant. Of the trinucleotide repeats (of at least four repeat units) AAT and AAC were the most frequent, followed by TTC. Excluding the dinucleotide repeats, the remaining 70 415 SSR loci were screened for the presence of sufficient forward and reverse flanking sequences suitable to design primers. This yielded 29 886 (42 %) PALs. Further filtering of these PALs by applying the most stringent criteria aimed at selecting single-copy loci yielded 4071 potentially amplifiable SSR loci.

Figure 1.

The 20 most frequent SSR motifs obtained, sorted according to frequency.

Polymorphism and amplification of SSR loci

A total of 136 SSR loci (117 randomly picked and 19 loci predicted to be potentially polymorphic as they appeared to have two different alleles in the sequence reads) were tested for amplification and degree of polymorphism in 10 randomly chosen individuals from different populations. Of the 117 randomly picked loci, 82 primer pairs amplified a high-quality PCR product, of which 37 (45 %) were polymorphic with a banding pattern that could be scored clearly (Table 1). Of the 19 primer pairs predicted to be polymorphic, 13 amplified bands of which 9 loci (69 %) were polymorphic, indicating a significantly higher rate of polymorphism (χ²test, P< 0.005) compared with randomly picked loci. The final set of 46 markers included 30 trinucleotide repeat markers, 4 tetranucleotide repeats, 6 pentanucleotides and 6 hexanucleotide repeats. The number of alleles across the polymorphic loci varied between 2 and 12 with an average value of 4.8 alleles in 10 genotypes. Several of the polymorphic markers with 10–12 alleles were TTC repeats. The heterozygosity per locus ranged widely from 0.10 to 0.89 (average 0.43). It is possible that, when used in larger populations, these markers will show higher estimates of Ho, and additional alleles may be found.

As shown in Table 1, most of the SSRs successfully amplified in B. pirrotae and B. popoviana (in the latter species the amount of DNA was insufficient to test all markers). Amplification, even if it is in the same size range as the alleles in B. papyrifera, is not proof that the marker is polymorphic, but heterozygosity (two different alleles in the expected range) is. Based on that criterion at least 19 of the 46 markers are polymorphic in B. pirrotae and at least 8 of 33 tested are polymorphic in B. popoviana.

Sequence assembly and annotation

The Illumina reads are the first genomic resource generated in the genus Boswellia. The repeat fraction was assembled based on k-mer frequency. This produced 49 576 contigs of repeats that were present at least 50× (median length 139 bp, mean length 224 bp, N50 238 bp, maximum length 21 153 bp, total sum = 574 Mbp). Next, 1533 contigs had blastn hits with RepBase, mostly with Copia (639 hits) and Gypsy (523) retrotransposons, alongside EnSpm (114), hAT (72), Satellite (29), TY (23), Harbinger (16), YPrime (14), Helitron (12) and SCTRANSP (3). Intermixed with these elements were hits to the ribosomal RNAs (LSU 56 hits, SSU 41) and also to Caulimoviridae viruses (11).

Using all data in a de novo assembly with SOAPdenovo, 444 927 contigs were obtained with a median of 375 bp, a mean contig length of 690 bp, an N50 of 1085 bp and a maximum contig size of 19 236 bp (total sum = 307 Mb genomic DNA sequence). The contigs >1000 bp were blasted against Genbank, and 65 467 were annotated with GO terms (Fig. 2; note that these are overlapping classes).

Figure 2.

Representation of ontology assignments of the B. papyrifera contigs. (A) The 31 086 GO terms of cellular components, (B) the 42 423 GO terms of molecular function and (C) the 54 256 GO terms of biological processes. Note that these are overlapping classes.

Terpene biosynthesis genes

Assefa et al. (2012) conducted a biophysical and chemical study on resins of Boswellia species with special emphasis on B. papyrifera. Using the list of identified components, eight contigs of the assembly were found, which represent part of genes of the terpene synthesis pathways, namely pinene synthase, limonene synthase (2×), isoprene synthase (4×) and gamma-terpinene synthase.

We also searched for the enzymes that are involved in terpenoid backbone biosynthesis (according to the Kyoto Encyclopedia of Genes and Genomes pathway database). Table 2 lists the enzymes of the mevalonate and non-metavolate (MEP/DOXP) pathways, the two pathways for the synthesis of terpenoid building blocks in plants, which were found among the annotation results. Two of the key enzymes of the MEP pathway, 2-C-methyl-d-erythritol-4-phosphate cytidylyltransferase (EC 2.7.7.60) and 4-(cytidine 5′-diphospho)-2-C-methyl-d-erythritol kinase (EC: 2.7.1.148), were not recognized in the set of scaffolds, but reciprocal tBlastx (at 1e-5) against these enzymes identified in Arabidopsis did reveal hits with, respectively, 3 and 2 contigs.