Genome and transcriptome sequencing characterises the gene space of Macadamia integrifolia (Proteaceae)

Genome sequencing and assembly

Approximately 37% of the assembled genome is identified as repetitive. As reported in most other plant species, long terminal repeats (Gypsy and Copia LTR) comprised the largest group accounting for approximately 29% of known repetitive elements and ~11% of the assembled genome (Fig. 1). Short and long interspersed repeats (SINEs and LINEs) accounted for ~18% while the majority of the identified repeats (~41.5%) were unclassified, lacking similarity to known repeats. In total, 98,114 perfect simple sequence repeat (SSR) motifs with di-, tri-, tetra-, penta- and hexanucleotide repeats were detected. Of these, 56,817 (57.9%) were dinucleotide repeats and consistent with reports for other plant species, the majority of these (58%) were AG/CT repeats [30]. In addition, there were 21,912 tri-, 11,262 tetra-, 5,045 penta- and 3,078 hexanucleotide repeats.

Transcriptome assembly, gene prediction and functional annotation of proteins

Transcriptome assembly using the quality controlled reads from three cDNA libraries (flower bud, young leaf and shoot) in Trinity de novo generated 298,030 contigs (transcripts) including different isoforms per contig. The transcripts had an N50 size of 1339 bp, mean transcript length of 823 bp, maximum transcript length of 17,814 bp and minimum transcript length of 224 bp (Table 1). Initial transcripts were clustered to generate a final set of 244,925 transcripts, which were used as one source of evidence in the evidence-based gene model prediction pipeline. Final annotation using MAKER pipeline and assembled transcripts produced 35,337 high-confidence gene models. Of these, 90.3% (31,908) were supported by expression values, FPKM (fragments per kilobase of transcript per million mapped reads) of 1 or more, and 87.6% (30,940) were supported by at least two RNA-seq reads. Although 3430 gene models lacked RNA-seq read support it is important to note that RNA-seq data was collected from flower, leaf and shoot tissue only. Over 78 and 74% of predicted proteins had at least one significant BLASTP hit (1E-05) against Nelumbo nucifera or Arabidopsis thaliana proteins respectively.

Core eukaryotic genes (CEGs) are 248 highly conserved genes understood to be present in virtually all eukaryotes in a reduced number of paralogs [31]. Among flowering plants, 959 single copy genes have been identified that are shared between Arabidopsis, Oryza, Populus and Vitis [32]. More than 84% of these single copy genes (809 genes) and 96% of CEGs (237 genes) had a significant BLASTP hit (1E-05) against the predicted macadamia genes. Assessment of annotation completeness with BUSCO (benchmarking universal single-copy orthologs) [33] indicates that the macadamia gene space contains 77.4% of the expected content. Using a 429 single-copy eukaryote gene set, 192 complete single-copy, 90 complete duplicated, 140 fragmented, and 97 missing universal single-copy genes were identified. This compares to 94.6% (23 missing) and 89.7% (44 missing) of the expected content in the high-quality genome assemblies of Eucalyptus grandis and Nelumbo nucifera respectively. In total, 19,794 macadamia genes were assigned to 33,291 InterProScan (IPR) domains and 39,925 GO terms. Predicted macadamia genes with a significant BLAST hit in KASS (KEGG Automatic Annotation Server) were assigned to 349 known metabolic or signalling pathways. The metabolic pathway (ko01100) contained the largest number of genes (826), followed by biosynthesis of secondary metabolites (ko01110, 386 genes), biosynthesis of antibiotics (ko01130, 188 genes) and microbial metabolism in diverse environment (ko01120, 147 genes).

Gene family analysis

Comparative genome wide analysis of orthologous genes was performed with OrthoVenn [34] to compare putative Macadamia integrifolia protein sequences with those of five other eudicot species including the core eudicots Arabidopsis thaliana, Eucalyptus grandis, Populus trichocarpa, Vitis vinifera and the basal eudicot Nelumbo nucifera. In total, 207,057 sequences from the six species were grouped into 23,778 clusters. Of these, 17,314 clusters contained at least two species and 1412 were single copy clusters containing one gene for each of the six species. There were 8743 orthologous gene clusters shared across all six species indicating their conservation within eudicots, while 1005 clusters containing 3168 genes were specific to Macadamia (Fig. 2). Macadamia and Nelumbo shared 587 gene clusters, the highest between any two species compared, consistent with their positions among basal eudicot families.

There was also evidence for an expansion of genes involved in terpenoid biosynthesis. In total, 78 macadamia gene models were functionally classified using Interpro as belonging to the terpene synthase gene family. Of these, 30 had high protein sequence similarity (1E-025) in BLASTP comparisons to Arabidopsis thaliana TPS-b monoterpene synthases. Among macadamia-specific clusters, significantly enriched GO terms included 25 predicted genes in six clusters involved in terpenoid biosynthetic process (GO:0016114, P = 0.036), and in particular biosynthesis of monoterpenes through geranyl diphosphate metabolic process (GO:0033383, P = 0.004) and monoterpene biosynthetic process (GO:0043693, P = 0.013). Monoterpenes, or C₁₀ isoprenoids are components of essential oils and fragrance in aromatic plants with roles in pollinator attraction, plant-plant interaction and defense with potential as pesticides and antimicrobial agents. While the functionality of these putative genes is yet to be tested, these results suggest that there may have been a lineage-specific expansion in macadamia of gene families involved in monoterpene synthesis.

Identification of candidate genes potentially involved in cyanogenic glycoside biosynthesis

In Macadamia, the cyanogenic glycoside (CG) dhurrin and its diglucoside derivative proteacin have been identified [23]. The metabolic pathways for cyanogenesis are best understood in Sorghum bicolor, Trifolium repens and Prunus spp. with three genes (CYP79, CYP71 and UGT85) encoding enzymes in the CG biosynthesis pathway from amino acid. Synthesis from specific amino acids is catalysed by cytochrome P450s and UDP-glucosyltransferase. Cyanogenesis occurs upon tissue disruption with catabolism involving a β-glucosidase and release of hydrogen cyanide (HCN) that is either catalyzed by a α-hydroxynitrile lyase (HNL) or occurs spontaneously at high pH (Fig. 3).

Evidence for expansion of plant defense-related gene families

Rainforests are among the oldest and most diverse ecosystems [39]. Australian subtropical rainforests in particular, are ancient refugia with high levels of plant species richness, endemism and rainfall [19]. Recent evidence suggests that insects and pathogens are instrumental in the maintenance of plant species diversity in rainforests [40]. Likewise, elevated predator-pathogen pressure is hypothesised to increase and diversify plant chemical defense systems. Plants have developed a wide range of defense systems to respond to the biotic stresses exerted by the predators and pathogens with which they have co-evolved [41, 42]. Expansion of the receptor-like kinase genes in particular is purportedly in response to fast-evolving pathogens [43]. Comparative genomic analyses suggests that there has been a lineage specific expansion in macadamia of gene families with similarity to Arabidopsis LRR receptor-like serine threonine-protein kinases EFR and FLS2. These encode proteins that play a key role in pathogen recognition and the activation of plant defense response [44, 45] and it has been demonstrated that Arabidopsis EFR enhances bacterium resistance in dicot and monocot transgenic plants including rice [46]. Further research is needed to identify the complete suite of macadamia plant resistance and defense genes and to determine whether polymorphism at sites on candidate genes is associated with resistance to co-evolving pathogens in macadamia as has been previously reported in Arabidopsis [47, 48]. Future growth in macadamia global production is expected following rapid expansion of cultivation and demand particularly in Asia. Germplasm collections, including clones of wild and domesticated trees have been established. These resources, along with wild populations undoubtedly contain genetic variants of interest for breeding, including improved yield, nutritional benefits, pest resistance and capacity to grow under variable climatic conditions. Insect herbivores and microbial pathogens are a major cause of yield reduction in macadamia production and identification of natural resistance would be of benefit for crop improvement.

Genes involved in cyanogenesis

Cyanogenesis is a plant chemical defense response to generalist herbivores involving the release of hydrogen cyanide following tissue disruption and hydrolysis of cyanogenic glycosides (CGs) [49, 50]. Endogenous recycling without cyanide release suggests that CGs serve additional biological roles including nitrogen and carbon supply at specific plant developmental stages [51] and there is evidence that intermediate compounds produced during biosynthesis of CGs have anti-microbial activity [52–54]. While relatively few plant species are cyanogenic, they are over-represented among food plants [25] and are common in the Proteacaeae, particularly in the subfamily Grevilleiodeae to which macadamia belongs [23, 24].

In macadamia, cyanide has been detected in seed, root, cotyledon and leaf tissue [22]. While levels in mature kernels of the commercial species M. integrifolia and M. tetraphylla are extremely low, they are much higher in the bitter mature kernels of M. ternifolia and M. jansenii. In almond Prunus amygadalus, bitterness of the kernel is determined by the content of the cyanogenic diglucoside amygdalin [55]. Intraspecific and temporal variation in cyanogenic capacity, and acyanogenic individuals have been reported in a number of cyanogenic plant taxa (e.g. [56, 57]). In white clover Trifolium repens, inheritance follows a Mendelian two-locus model. The Ac/ac (CYP79D) gene controls production of cyanogenic glycosides, and the Li/li (cyanogenic β-glucosidase) gene controls their hydrolysis [58]. There is an apparent selective advantage for acyanogenic individuals in colder climates, and polymorphism is maintained within populations through recurrent gene deletions over time [59].

We identified macadamia homologues with high sequence similarity to genes encoding enzymes involved in CG biosynthesis in other cyanogenic plants including Sorghum bicolor and Trifolium repens. Based on the relatively high RNA-seq expression values in green tissue, six homologues to genes encoding the enzymes CYP79, CYP71, UGT85 and β-glucosidase are probable candidates in macadamia to target for further analysis. The discovery of candidate cyanogenesis genes in macadamia is likely to be an important step in facilitating the utilization of the smaller tree species M. ternifolia and M. jansenii into breeding programs to reduce tree size while retaining kernel edibility. In previous studies, 28% of the Proteaceae species tested were cyanogenic. This compares to 4.5% of 401 species from 87 families in Australian rainforests, and 4% of Eucalyptus species [24, 60]. The high proportion of cyanogenic plants in Proteaceae indicates that cyanogenesis is an important defense strategy in this family. Further work is planned to validate candidate genes, screen wild macadamia germplasm for natural variants and investigate the interaction between pest resistance, climatic variation and cyanogenesis in macadamia and more broadly across the Proteaceae.

Plant materials

Fresh plant tissue was collected from a Macadamia integrifolia, cultivar 741 ‘Mauka’ individual from the Macadamia Varietal Trial plantation M2 at Clunes, New South Wales, Australia and stored at -80 °C. A voucher specimen is deposited in the Southern Cross University herbarium [accession PHARM-13-0813]. Prior to DNA and RNA extraction, leaf tissue was frozen in liquid nitrogen and ground using a tissue lyser (MM200, Retsch, Haan, Germany).

Genomic DNA isolation and sequencing

Total genomic DNA was extracted using a DNeasy Plant Maxi kit (Qiagen Inc., Valencia, USA) for all DNA sequencing with the exception of mate pair (MP) library sequencing where DNA was extracted using a CTAB-based method developed for next-generation sequencing [61]. DNA was quantified using a Qubit dsDNA BR assay (Life Technologies, Carlsbad, USA). Genomic DNA was sheared using a Covaris S220 focused-ultrasonication device (Covaris Inc., Woburn USA). Paired-end libraries (PE) with average insert sizes of 480 and 700 bp and an 8 kb MP library were prepared using Illumina TruSeq DNA Sample Preparation kit v2 following manufacturer’s instructions (Illumina, San Diego, USA). Fragment size distribution and concentration were determined using a DNA 1000 chip on a Bioanalyser 2100 instrument (Agilent Technologies, Santa Clara, USA). PE and MP libraries were sequenced with Illumina GA IIx (150 x 2 cycles) and HiSeq 2000 (100 x 2 cycles) instruments respectively.

Genome assembly and scaffolding

Paired-end sequence reads were trimmed to remove low quality bases and adapter sequences and de novo assembled using CLC Genomics Workbench (CLC) version 6.5 (CLC Bio, Aarhus, Denmark) that has been used in the assembly of plant genomes including Norway spruce Picea abies [62] and barley Hordeum vulgare [63]. CLC de novo assembler, which utilizes de Bruijn graphs, was used for assembly of Illumina PE reads with the option to map reads back to contigs following previously described parameters [4]. MP reads were also trimmed to remove low quality bases and adapter sequences. We observed very high proportion (>90%) of duplicated MP reads, presumably PCR duplicates, which were filtered using CLC. Genome assembly was performed in the following two steps: preliminary contig assembly using PE reads in CLC, followed by assembly of sequence contigs and filtered high quality MP reads using the scaffolding program SSPACE to obtain a final set of scaffolds [64]. Genome size was estimated based on k-mer analysis and depth of sequencing [29].

Repetitive sequence analysis

RepeatModeler and RepeatMasker programs were used to identify repeats [65]. Putative repetitive sequences were identified using the RepeatModeler program with default parameters. In parallel, known repetitive sequences were identified using the RepeatMasker program with the latest release of RepBase curated repeat libraries [66]. Searches for simple sequence repeats (SSRs) were conducted using SciRoko [67] software with default parameters and ‘MISA’ mode.

RNA extraction and transcriptome sequencing

To enable assembly of the transcriptome of macadamia, three tissues (leaf, shoot and flower) of cultivar 741 ‘Mauka’ were selected for deep RNA sequencing (RNA-seq). Total RNA was isolated from frozen tissue using Ambion Plant RNA Isolation Aid prior to extraction using an Ambion RNAqueous Kit following manufacturer’s recommendations (ThermoFisher Scientific, Waltham, USA). Libraries were prepared with Illumina TruSeq Stranded mRNA Library Preparation Kit and PE sequenced with an Illumina HiSeq 2500 (100 x 2 cycles).

Transcriptome assembly

Quality control of tissue specific transcriptomic reads involved removal of low quality sequences, adapter sequences and empty reads using BBMap tools (sourceforge.net/projects/bbmap/). Retained high quality clean reads were assembled using the Trinity de novo transcriptome assembly program (version 2.0.2) with default parameters [68]. The Trinity de novo assembly pipeline consists of three different modules Inchworm, Chrysalis and Butterfly. Inchworm assembles short reads into unique sequences of transcripts. Chrysalis clusters the Inchworm transcripts and constructs de Bruijn graphs for each cluster where each cluster represents the full transcriptional complexity for a given gene. Butterfly then processes the individual graphs in parallel, tracing the paths that reads and pairs of reads take within the graph, ultimately reporting full-length transcripts for alternatively spliced isoforms, and resolving transcripts that corresponds to paralogous genes. The initial transcripts were clustered using the CD-hit-est [69] to generate final set of transcripts, which were used as one source of evidence in the evidence-based gene model prediction pipeline.

Gene prediction and annotation

Annotation of gene models was conducted using MAKER (version 2.31.8) which is an evidence-based gene model prediction pipeline [70]. MAKER combines the power of protein and Expressed Sequence Tag (EST) based homology with ab initio gene predictions to produce polished gene annotations. Trinity assembled transcripts were used as reference ESTs, and proteins of Nelumbo nucifera and Arabidopsis thaliana were used as reference proteins [35]. Macadamia scaffolds were first repeat masked using RepeatMasker [65]. To obtain the homology based genes MAKER aligned reference ESTs and proteins using Blastx [71] and exonerate [72] against the macadamia scaffolds. Ab initio gene predictions were made by Augustus [73] and SNAP [74] gene prediction programs. MAKER created the final gene set by combining the evidence based and ab initio predictions.

Functional annotation of proteins

Predicted protein coding genes were functionally annotated based on protein signatures and orthology relationships. Similarity search was performed against release (03-2015) of UniProt Swiss-Prot proteins. Functional domains, gene ontology (GO) terms, GO accessions were searched against InterPro using InterProScan software [75]. Functional and gene ontology (GO) domains were assigned using InterProScan as described in [76] with default parameters. InterProScan integrates a collection of protein signature databases including BlastProDom, HMMPfam, HMMSmart, HMMTigr, ProfileScan, HAMAP, PatternScan, SuperFamily, TMHMM, HMMPanther, Gene3D and Phobius. To inform biological interpretation of macadamia gene function, KEGG (Kyoto Encyclopedia of Genes and Genomes) reference pathway database was used to map macadamia genes to defined pathways [77]. The KASS (KEGG Automatic Annotation Server) was used to assign genes to metabolic pathways using BLASTX with an E-value cutoff of 1E-05 [78]. Tests for annotation completeness were conducted using BUSCO [33] with the eukaryote 429 gene set and compared results to those of the Eucalyptus grandis and Nelumbo nucifera genomes.

Comparative genomic analysis and gene family identification

Protein sets of five plant species including core eudicots Arabidopsis thaliana, Eucalyptus grandis, Populus trichocarpa, Vitis vinifera and basal eudicot Nelumbo nucifera, were downloaded from respective public repositories. Along with the predicted macadamia proteins they were uploaded into the OrthoVenn web server for identification and comparison of orthologous clusters [34]. To identify orthologous groups OrthoVenn employs the OrthoMCL Markov clustering algorithm, although unlike OrthoMCL it employs UBLAST for the all-against-all similarity search, which is ~350 times faster than conventional BLAST [79]. Following clustering, orthAgogue [80] is used for the identification of putative orthology and inparalogy relations. To deduce the putative function of each ortholog, the first protein sequence from each cluster is searched against the non-redundant protein database UniProt [80] using BLASTP. Pairwise sequence similarities were determined among protein sequences of all species with a BLASTP E-value cut-off of 1E-05 and an inflation value of 1.5 for MCL. To test the quality and completeness of the gene space assembly of macadamia we identified orthologous clusters from analyses in OrthoVenn with Swiss-Prot hits to proteins reportedly involved in CG biosynthesis and activation, and conducted BLASTP searches of the macadamia candidates. In addition, reciprocal searches of protein sequences for five enzymes (CYP79, CYP71, UGT85, β-glucosidase and HNL) involved in CG biosynthesis from known cyanogenic plants were conducted with a DeCypher Tera-BLASTP search against all macadamia gene models.