A manually annotated Actinidia chinensis var. chinensis (kiwifruit) genome highlights the challenges associated with draft genomes and gene prediction in plants

Assembly of a second genome of Actinidia chinensis var. chinensis

To generate a new Actinidia chinensis genome, a diploid F3 sibcross individual Red5, with a predicted inbreeding coefficient of 37.5%, was chosen for sequencing (Fig. 1). An anytag-based assembly of paired-end Illumina reads generated 46,117,212 fragments with an N50 of 275 bases [26]. Assembly of these fragments using a long insert library (Roche 454 GS-FLX – 4 kb) using Newbler produced 39,868 contigs. Subsequent stepwise scaffolding using SSPACE2 [27] with Illumina long-range insert libraries of 4, 9 and 13 kb yielded 39,825, 8688 and 3887 scaffolds, respectively. After two iterations of gap closure the final assembly consisted of 3887 scaffolds with a total length of 550.5 Mb. The N50 was 623.8 kb with L50 of 240 scaffolds, an N90 of 140.7 kb with L90 of 941 scaffolds, and 3.57% N content, with the longest scaffold being 4.43 Mb.

The genetic linkage map of Scaglione and colleagues [15], which included markers from Fraser and colleagues [28], augmented by BLAST walking comparison between Red5 and scaffolds of ‘Hongyang’ enabled anchoring of 2727 scaffolds (Table 1) comprising 547.9 Mb, to 29 linkage groups. The remaining 1206 unanchored scaffolds containing 5.91 Mb with an N50 of 5.36 kb were concatenated to form a composite entity hereafter referred to as ‘Chr30’ for the purposes of subsequent manual annotation of the entire genome sequence. Note, upon submission to NCBI Genbank scaffolds assigned to ‘Chr30’ were submitted as individual (non-concatenated) scaffolds according to NCBI Genbank submission policy. Estimates of genome size based on K-mer analysis indicated a genome size of 705 Mb (preQC) or 742 Mb (jellyfish). These align with estimates from flow cytometry [29] that report the genome to be 758 Mb in size. The assembly therefore represents ~ 73% of the estimated genome size with 98.9% of assembled scaffolds (72% of the estimated genome size) assigned to chromosomes. This is a considerable improvement from the original ‘Hongyang’ draft genome, which had 164 Mb unassigned (Fig. 2).

Chromosome	Size (Mb)	Number of Scaffolds	Scaffold N50 (Kb)	Longest Scaffold (Mb)	Manually Annotated Genes
					Number of Genes	Number of Transcripts	Genes Per 100 Kb
1	18.6	94	595.4	1.51	1125	1133	6.06
2	14.6	83	536.6	1.41	1091	1094	7.46
3	21.7	92	1154.9	2.96	1726	1729	7.94
4	13.8	75	420.4	2.39	781	781	5.67
5	18.6	113	655.2	4.44	960	976	5.16
6	17.4	109	474.8	1.56	1059	1069	6.09
7	20.0	99	666.0	1.82	1068	1075	5.33
8	26.1	134	779.6	1.68	1448	1459	5.55
9	16.6	66	598.4	1.95	1118	1123	6.74
10	19.3	111	413.0	0.95	986	994	5.10
11	16.9	74	520.3	1.26	1119	1137	6.63
12	19.2	122	480.2	1.42	1053	1054	5.49
13	19.5	77	926.8	1.65	1388	1402	7.11
14	17.9	75	651.2	2.32	1106	1129	6.18
15	15.9	67	707.9	2.13	1106	1121	6.94
16	23.8	152	388.5	1.46	1252	1255	5.26
17	17.4	100	387.2	0.91	931	934	5.34
18	20.7	85	850.4	2.80	1213	1219	5.85
19	15.4	117	345.9	2.12	653	656	4.24
20	17.9	89	512.5	1.53	1055	1056	5.88
21	17.3	71	871.2	2.16	1046	1064	6.04
22	18.9	105	462.5	1.42	1092	1096	5.76
23	27.7	79	802.6	2.15	2324	2327	8.39
24	17.8	59	1032.0	3.02	1198	1201	6.72
25	19.6	83	966.1	2.45	1008	1010	5.14
26	20.4	79	1086.4	3.41	1237	1247	6.07
27	21.0	112	461.8	1.27	1013	1040	4.83
28	15.8	75	964.5	2.80	1011	1012	6.38
29	18.0	130	402.3	0.69	993	1002	5.53
Total for Chrs 1–29	548.0	2727	575.6	4.44	33,160	33,395	6.05
Unassigned 30	5.9	1206	5.4	0.66	97	97	1.64
Total	553.9	3933			33,257	33,492	6.00

Evaluation of genome assembly

To evaluate genome assembly accuracy we assessed mapping back of paired end reads to the assembly and compared the assembled contigs for 22 clones from a BAC library of A. chinensis Red Female 1 (Fig. 1). These contigs resulted from sequencing using a different technology (454) and different assembly path (Newbler). Rates of discordant alignment of input paired end reads mapping once to the whole genome sequence was low (0.33 to 2.24%) (Additional file 1). Alignment of these BAC clone contig assemblies to chromosome 25 supported the assembly of Red5. The alignments (Additional file 2) show a close correspondence between the BAC clone assemblies and the whole genome assembly of Red5 in this region of chromosome 25. Additionally, the alignment of read pairs from the 9Kb LIMP library was assessed and visualised using hagfish_blockplot from the ‘hagfish’ software (https://github.com/mfiers/hagfish/). The majority of alignments displayed green indicating the read pairs from the 9Kb LIMP library aligned to the chromosome sequences with the default bounds determined by ‘hagfish’ (Additional file 3). As expected for assembly from short read data there were also regions depicted in pinkish-red suggesting that the mate pairs aligned to the genome in these regions outside the expected distances. Such regions will occur for example when repeats are compressed into a consensus leading to a compression in the whole genome assembly sequence relative to the physical genome sequences.

Global analysis comparing the new genome with itself revealed areas of similarity among different regions of the chromosomes (Fig. 3a). When these duplicated chromosomes were examined more closely it was found that many appeared to have Robertsonian–like centromeric translocations. This could be clearly seen, for example, in the duplicated Chr1, which showed homology to half of Chr8 and half of Chr9 (Fig. 3b). The paired chromosomes were arranged sequentially (Fig. 3c) and all but two had at least one translocation event. The only non-translocated chromosomes were Chr4 (homeologous to Chr21) and Chr2 (homeologous to Chr3) (Fig. 3c).

Manual curation of the predicted gene models

To develop gene models for the new genome, a WebApollo tool [21] was populated with the new genome sequence and the following tracks of evidence were added:

1. 1)

   A scaffold assembly quality track which identified single N insertions, indicative of construction anomalies, and N repeats that marked the edges of scaffold boundaries

    
2. 2)

   A repeat masker track that identified repetitive elements, including transposons

    
3. 3)

   Computational gene prediction tracks consisting of the original ‘Hongyang’ gene models, and an ab initio gene prediction scan of the Red5 genome

    
4. 4)

   ESTs from published EST sequencing libraries, and 12 tracks of RNA-Seq from diverse tissues from Red5 (Table 2).

    

An international consortium of annotators synthesised this information to produce a new gene model for each gene along the genome.

During the annotation process, some models were relatively simple to predict while others were more difficult. Easy-to-annotate models had strong RNA-Seq support, clear intron-exon structures and a gradation in the number of RNA-Seq reads at the beginning and the end of the gene. These genes also often had an accurate computationally predicted gene model and had a good BLASTP match from GenBank that covered the majority of the gene. Conversely, there were many regions of the genome that were very difficult to annotate. These harder-to-annotate regions had conflicting, patchy, or even no RNA-Seq evidence, combined with conflicting or absent computational gene models (Fig. 4). The complexity of these gene regions was often confounded by genome structure caused by repetitive elements and/or anomalies in the genome sequence construction, as observed by a single “N” or scaffold boundaries in the quality track. Within the genome, the whole spectrum of combinations of these challenges were identified with occasional loss of open reading frame caused by the anomalies. To address this variation in confidence for each model, a quality tag was added to each of the gene models. A strongly supported high quality gene model was given a Q2 tag, while a model that was suspected to be incomplete or contain other errors was given a Q1 tag. Further tags based on whether the BLASTP alignment suggested it was full length (F) or a partial (P) gene were also added. A Q0 tag was given to computationally predicted genes that exhibited peptide homology to a GenBank protein, but no RNA-Seq evidence. Following the first pass of manual annotation, it was found that there was a high degree of variability among the 93 annotators in their interpretation of the evidence to create the models, especially for loci where evidence was ambiguous, sparse, or apparently contradictory. To address this variation, a second pass of the genome was undertaken by a smaller group of ‘expert’ annotators to standardise the annotation.

The manual annotation identified 33,044 gene loci, with 33,123 protein isoforms. The mean coding length of annotated genes was 1278 bases, while the mean locus length was 6071 bases (range 282–92,048 bases). Of the 33,123 loci, 32,967 (99.5%) were manually annotated with a single isoform, while 76 were annotated with two isoforms. One locus was annotated with four isoforms respectively. It should be noted that for loci with potentially more than one isoform individual isoforms were only annotated where they were cleanly defined from the RNA-Seq evidence. Otherwise a single model was submitted by annotators. Therefore it is likely that there will be further isoforms defined as further evidence is obtained. Of the 33,123 loci, 2485 were annotated by the community annotators as partials, with 1057 annotated as having 5′ truncations, and 817 as having 3′ truncations, while 611 loci were annotated as being partial without an indication of truncation direction. The total number of coding sequence regions (CDS) in annotated genes was 181,135 and contained 42.3 Mb of genome sequence (5.58% of the estimated genome size). There were 6514 loci (19.7%) annotated as containing a single CDS. The mean length of CDS regions in single CDS loci was 1004 bases. In comparison, the mean CDS length for all loci was 233 bases. The mean intron length within the coding regions of multiple exon-containing loci was 886 bases. The minimum intron length was 22 bases. The corresponding quality scores were assessed and 83.16% of gene models had a Q2 score, 14.62% had a Q1 score and 2.18% had a Q0 score.

Since genomes are always evolving, it is likely that different genomic structures will be observed across Actinidia species. Additionally, new improved versions of the genome annotation will be developed that incorporate the last 6 Mb of unassigned fragments. To address this, we named our genes each with a unique name that will be enduring and independent of chromosome location. We have used the nomenclature AccXXXXX with splice variants appended as ‘.1’ and ‘.2’ etc., and then appended chromosome location as a descriptor which can subsequently be changed, without changing the unique name of the gene. To assist the kiwifruit research community we provide a conversion table for the best reciprocal matches between our gene set and that within the Kiwifruit Information Resource (KIR) [18] as Additional file 4 and will also make this data available via our Git Repository (https://github.com/PlantandFoodResearch/Red5_WGS_Manual_Annotation).

In addition to the location, the sequence length, sequence type (cds, cdna, peptide), and manual annotation quality score were appended along with an internal database identifier and functional description. The manual quality score also contained the F and P labels based on BLAST match, so a point to note, is a gene with a manual quality score of 2 score (good RNA-Seq support) but a BLAST alignment indicating truncation could be scored as 2P with the direction of truncation indicated by using a 5 or 3 suffix, for example a 5 prime truncation would be scored ‘2P5’.

Comparison of the Red5 gene set with the ‘Hongyang’ gene annotation

The manually annotated gene set was compared with the 39,040 published ‘Hongyang’ gene models originally published [14] (hereafter termed original ‘Hongyang’ model set) as well as to the 39,761 revised gene annotations [18] (hereafter termed revised ‘Hongyang’ model set). As ‘Hongyang’ is a different cultivar, polymorphisms are expected. To get a more accurate comparison, predicted protein sequences were used. When Red5 was used as query against the original ‘Hongyang’ model set as a database it was found that only 1973 (~ 6%) of the protein sequences for the 33,123 isoforms were identical in sequence and length to a ‘Hongyang’ predicted protein model [14]. We also detected instances where a Red5 model was perfectly contained within a ‘Hongyang’ model suggesting either the Red5 model is truncated, the ‘Hongyang’ model is over predicted or there is a genotypic difference between the two genotypes. The reverse situation where a ‘Hongyang’ model was perfectly contained in a Red5 model was also encountered. 882 (2.67%) of Red5 proteins were perfectly contained within a longer sequence of an original ‘Hongyang’ model while 828 (2.51%) Red5 proteins perfectly contained the sequence of an original ‘Hongyang’ model. When the revised ‘Hongyang’ model set [18] was employed as the database 3114 (9.4%) Red5 protein sequences were found to be identical in sequence and length to a protein within the revised ‘Hongyang’ model set while 927 and 1007 Red5 proteins respectively either were perfectly encapsulated within a revised ‘Hongyang’ model or perfectly encapsulated a revised ‘Hongyang’ model. We repeated the analysis in the reverse direction. When using Red5 as the database, 42% of the original ‘Hongyang’ model proteins and 54.2% of the revised ‘Hongyang’ model proteins possessed a match with identity of 90% or greater, showing a considerable number of genes have been changed firstly in the revised annotation and secondly in the manual annotation process.

Comparison of the Red5 to original ‘Hongyang’ models identified 1958 original ‘Hongyang’ models had identical sequence and identical length to the corresponding Red5 model. A further 1261 original ‘Hongyang’ models possessed identical protein sequence to a Red5 model but were shorter in length than the Red5 model while for 576 original ‘Hongyang’ models the reverse was true. As expected 3114 of the revised models were found to have identical sequence and identical length to the corresponding Red5 model. A further 1685 revised ‘Hongyang’ models possessed identical protein sequence to a Red5 model but were shorter in length compared to the Red5 model while for 553 revised ‘Hongyang’ models the reverse was true. To examine the relationship with less than perfect matching best reciprocal BLASTP matches between Red5 and ‘Hongyang’ protein datasets identified for 19,179 proteins [14] and 21,479 proteins [18]. When the lengths of predicted proteins identified as best reciprocal BLASTp matches were compared, 5542 and 4700 proteins respectively within original and revised ‘Hongyang’ genes, respectively, possessed a longer protein sequence length than the Red5 model. By comparison 13,482 and 16,551 Red5 protein were longer than their best reciprocal BLASTp match counterparts from original and revised ‘Hongyang’ model sets respectively.

Within both the original and revised ‘Hongyang’ gene sets 148 and 113 models were completely missing from the Red whole genome sequence. The identifiers for these models are listed in Additional file 5. For 1195 original and 587 revised ‘Hongyang’ models lacking a best reciprocal BLASTp match we found the CDS for these models to be encapsulated in the UTR regions of Red5 models. For a further 379 original and 82 revised ‘Hongyang’ models lacking protein:protein matches to Red5 proteins we found their CDS to overlap the 3 UTR of a Red5 model while for 362 original and 77 revised ‘Hongyang’ models the CDS was found to overlap the 5 UTR of a Red5 model. A further 3534 and 2034 models from the original and revised ‘Hongyang’ sets respectively were completely present in the whole genome sequence of Red5 but possessed no protein match to a Red5 gene model and did not align to a UTR region of a Red5 model. To identify if these models were missing from our annotation set due to lack of support from RNA-Seq evidence we merged BAM files for RNA-Seq libraries previous aligned to the Red5 whole genome sequence for purposes of assisting manual annotation. The CDS for the ‘Hongyang’ gene set of Yue and colleagues [18] was aligned to the Red5 whole genome sequence using GMAP [31] (version 2017–06-20), the resultant GFF3 output converted to Simplified Annotation Format (SAF) and RNA-Seq read counts to these features extracted using featureCounts [32].

Of the 2034 revised ‘Hongyang’ models that perfectly aligned to the Red5 whole genome sequence but for which there was not protein match of any kind, 535 aligned to regions of the Red5 genome where there was no aligned RNA-Seq and thus would have been unlikely to be manually annotated as a result. To further examine the RNA-Seq alignment of the remaining 1499 revised models, the base coverage on each chromosome was extracted using bedtools genomecov (v2.21.0) [33]. A perl script was then used to convert genomecov’s chromosome base by base coverage to a bitmap (0 for no coverage at base position, 1 for coverage) and an array of coverage values for each exon of each of the ‘Hongyang’ 1499 revised alignments models. These were filtered to identity models incomplete coverage and models with coverage across their entire match regions. Of the 1499 revised models examined 131 possessed exons which were not supported by Red5 RNA-Seq while a further 372 revised ‘Hongyang’ models possessed RNA-Seq coverage of less than 5 reads per base. Given that we used all Red5 RNA-Seq combined in this analysis while annotators examined evidence library by library it is possible that these were not annotated due to inconsistent coverage across evidence libraries. The remaining 1069 regions represented 967 individual revised models each with RNA-Seq coverage of greater than 5 reads per base. We have provided a list of those revised models in Additional file 6 including their locations in the Red5 whole genome sequence and the average number of RNA-Seq reads aligned on a per base in each CDS.

Using the new manually annotated gene models, two gene families were investigated further to assess whether new annotation had missed any genes. Given the reported poor annotation of EXPANSIN (EXP)-like genes [17], the genome sequence was translated into all six translation frames and used to identify regions that had homology to Arabidopsis EXP, and EXPANSIN-LIKE (EXPL) protein sequences. Fifty three chromosomal regions were identified, and 41 of these had a new manually annotated gene model. Of the 12 remaining gene models, six were partial regions of homology from which no gene models could be generated. Two regions each coded for a possible full length gene even though there was no associated computer-predicted model and no RNA-Seq evidence to support them. These two genes were scored Q0 and added to the gene list, producing a possible 47 EXP genes. Comparison of these 47 EXP genes with the ‘Hongyang’ gene models showed that six were identical to the ‘Hongyang’ gene models (Fig. 5a, Additional file 9), 18 were partially supported by an original ‘Hongyang’ gene model and 15 by a revised ‘Hongyang’ gene model. 29–32 EXP models were new and interestingly the majority of these were EXP genes with the EXPL more accurately predicted. A second gene family, the ACC SYNTHASE (ACS) genes, was also investigated [34]. In total 16 translated chromosomal regions showed homology to ACS proteins. Of these 16, 14 had a manually annotated gene model associated with them (Fig. 5b, Additional file 9). For the two regions that did not have a manually assigned gene, no convincing gene models could be derived. Two genes, ACS12 and ACS13, had duplicate gene models that were 100% identical to each of the models at the nucleotide level, with ACS12 and ACS12R on two separate scaffolds on the unassigned chromosome (Chr30), suggesting these are possibly allelic. The other pair (ACS13 and ACS13R) were found as a tandem duplication on Chr12. When the 14 ACS genes were compared with the original ‘Hongyang’ gene models, the 12 unique genes were all found; nine of these models were identical to the manual annotation, and three had a single exon/intron boundary difference, these numbers were unchanged in the revised ‘Hongyang’ models. The EXP study suggests that it is likely, even after manual annotation, that there are still unannotated genes of other gene families in the genome, especially in gene families of similar structure to the EXP genes i.e. smaller and computationally hard to predict. However, the new manually annotated genes do contain a more comprehensive list than the previous computationally generated gene lists.

At a chromosomal level, the gene density along the chromosomes varied, with less well populated regions associated with lower recombination rates (Fig. 6a). These lower density regions have been previously linked to centromeric regions [35, 36], and are also often associated with the translocation cut sites, consistent with the Robertsonian translocation. Especially clear are translocations of Chr1 to Chr9 and Chr8, where estimations of centromere locations can be made based on gene density (Fig. 6b). In other chromosomes such as Chr6 it is less clear, suggesting that within the same chromosome there are a number of translocation events that have occurred. Fig. 6c demonstrates that these other non-Robertsonian translocations are evident in Chr19, where low gene density regions flank a region of homology to Chr6 with a high gene density (Fig. 6c). As well as a WGD and subsequent translocations there are a considerable number of localised duplicated regions. Indeed, 1572 sites in the genome contained tandemly duplicated genes, representing 9.43% of annotated genes.

Plant material

Two F₁ diploid A. chinensis Planch. var. chinensis from an open pollinated red-fleshed fruiting mother were screened with 8 microsatellites (Additional file 11) to ascertain that they were true siblings. These were crossed and F₂ offspring, a female with red fruit (Red Female 2) and male were selected for further crossing. Forty F₃ progeny were sown, with each having an inbreeding coefficient of 0.375 (Fig. 1). A red fruiting female (Red5) was chosen for genome sequencing. For gene expression analysis (RNA-Seq), different tissues were harvested from mature Red5 plants to encapsulate a diversity of expression (Table 3).

DNA isolation and sequencing

Nuclear DNA was isolated from leaf tissue of Red5 using nuclei enrichment and DNA extraction as described by Naim and colleagues [45]. DNA was sheared to an insert size of either ~ 160 bp or ~ 240 bp and prepared for 100 base paired-end sequencing along with 100 base long-insert mate-paired-end (LIMP) libraries with average insert sizes of 4, 9 and 13 Kb and sequenced on Illumina HiSeq2000™ (Illumina Inc. San Diego, CA, USA) at the Australian Genome Research Facility (AGRF - Brisbane), according to the manufacturer’s instructions. A 4 Kb insert library was also prepared for paired-end sequencing by Life Science (Roche) 454 GS-FLX. Cyclically corrected sequences from a small number (6) of PacBIO SMRT cells (45 Mb per cell) were also included during gap closing.

RNA isolation and sequencing

RNA was extracted using the method described in Chang and colleagues [46]. RNA samples were quantified and sample purity was verified by using a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific). RNA integrity was checked by an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). RNA was supplied to Macrogen Inc. (Seoul, South Korea) for standard RNA-Seq preparation and sequenced using the Illumina HiSeq2000™ yielding either single or pair-end RNA-Seq reads.

Bacterial artificial chromosome (BAC) library sequencing

A total of 11,520 clones from a BAC library made from nuclear DNA of the F1 mother (Red Female 1) were selected for sequencing. The re-arrayed BAC clones were grown in 96-well plates containing 1.2 mL LB liquid medium with 12.5 μg/mL tetracycline at 37 °C in an orbital shaker at 180 rpm for 16 h. The bacterial cells were harvested at 3000 rpm at room temperature for 30 min in a benchtop centrifuge. The BAC DNA was extracted using a plate based alkaline lysis method [47] and dissolved in 150 μL 28 mM Tris-HCl pH 8, 1 mM EDTA, 0.6 mM cresol red (to provide a visual aid for robotic transfers). Each BAC was individually barcoded using an in-house method BACRB (details of which can be supplied on request to the corresponding author). The barcoding oligonucleotides (BioSearch Technologies, Novato, CA, USA) were dissolved in TE buffer (10 mM Tris-HCl pH 7.5, 1 mM EDTA) to a final concentration of 50 pmol/μL. Approximately 15 ng of BAC plasmid (5 μL) was randomly tagged by primer extension pre-amplification PCR (PEP-PCR) in 50-μL reactions using AccuPrime™ Taq DNA polymerase High Fidelity system as per the manufacturer’s instructions (Life Technologies Corporation, Carlsbad, CA, USA) and 10 pmol of a corresponding BACRB oligonucleotide. The DNA was denatured at 94 °C for 2 min and amplified using 50 cycles of: 94 °C for 40 s, a two-step annealing strategy (30 °C for 2 min, ramp at 0.1 °C/s, 48 °C for 4 min), 68 °C 60s, followed by a final extension at 68 °C for 7 min. The randomly tagged BAC amplicons of each clone were amplified with the tailed BACRB oligonucleotide in a second touchdown (TD) amplification process. One microlitre of randomly tagged BAC amplicon was added to 19 μL of PCR mix (0.6 M trehalose, 40 mM Tris-HCl pH 8, 20 mM KCl, 20 mM (NH₄)2SO₄, 10 μg BSA, 0.5 mM MgSO₄, 0.2 mM dNTP, 2 pM tailed-BACRB oligonucleotide, and 0.25 unit Platinum® Pfx DNA polymerase (Life Technologies)). The touchdown amplification was performed as follows: 94 °C for 2 min, 1 cycle; (94 °C for 30s, TD 60–50 °C for 30s, 68 °C for 1 min), 20 cycles; (94 °C for 30s, 50 °C for 30s, 68 °C for 1 min), 10 cycles; and final extension at 68 °C for 7 min. The barcoded samples from each 384-well plate were pooled, concentrated, and analysed by agarose gel electrophoresis. A barcoded TruSeq library was prepared for each plate pool (30 barcoded libraries). A super pool was prepared by combining 10 barcoded libraries of plate pools. A total of three super pools were obtained, and each one sequenced separate lanes on single end mode, at Macrogen Inc. (Seoul, South Korea). The three lanes generated 576.35 million reads, comprising 58.2 Gbp.

For Roche 454 GS-FLX sequencing, 50-mL cultures were grown and extracted as described above. The DNA pellet was dissolved in 50 μL buffer (10 mM Tris-HCl, 1 mM EDTA, pH 7.5) and sent to be sequenced by 454 GS-FLX at Macrogen Inc. (Seoul, South Korea).

Assembly

A “PseudoSanger”-like approach [26] was used to assemble two paired-end read libraries with stepwise decreasing insert size (240 and 160 bases, respectively). The libraries yielded 169,008,438 and 170,367,691 read pairs, respectively. Prior to assembly, reads were error-corrected using the error correction tool from the ALLPATHS-lg assembler [48] yielding 159,232,897 and 167,054,602 corrected read pairs, respectively. These reads were also used to estimate genome size using both preQC (https://github.com/jts/sga/wiki/preqc) from sga [49] and jellyfish (version 1.1.10) [50]. Error-corrected reads were assembled using anytag (version 2.5.2) [26]. Anytag yielded 46,117,212 fragments with a minimum length of 81 bases, maximum length of 450 bases and N50 of 275 bases. These fragments were assembled using Newbler 2.9 (Roche 454, Bradford, Connecticut, USA) with settings “-m –large –het –cpu 32”. Also included in the Newbler assembly were 1,209,245 paired end sequences from a 4-kb insert library sequenced using 454 GS-FLX pyrosequencing. UniVec_Core (ftp://ftp.ncbi.nlm.nih.gov/pub/UniVec/UniVec_Core) and common sequencing primers were used for vector trimming during Newbler assembly, while the sequence of Escherichia coli DH10B (NC_010473.1) was used for screening for bacterial contamination.

Reads from Illumina LIMP libraries were trimmed to 36 bases and redundant pairs removed prior to use using a custom perl script. The resulting read pairs were then used to scaffold the Newbler assembly further using SSPACE2 [27], which was then followed by two iterations of gap closure using GapCloser (v1.12; http://soap.genomics.org.cn/about.html). Two additional Illumina paired-end read libraries, not used in the anytag/Newbler-based contig assembly, were also employed during gap closure along with the assembly input read libraries. These two additional libraries possessed average insert sizes of ~ 240 bp and ~ 300 bp with read lengths of 75 bp (trimmed to 64 bases) and 150 bp (trimmed to 109 bases) respectively. Gap closure using Illumina-derived sequence reads yielded a reduction in ambiguities (N) within the assembly from 17.3 to 3.57%. Cyclically corrected sequences extracted from 6 SMRT™ cells (PacBIO) containing 35–45 Mb per cell were compared with the assembly using BLAT [51] and the resulting alignments used to further infill 199.01 Kb of genome assembly gaps. This assembly process yielded 3387 scaffolds containing 554 Mb with a minimum size of 1997 bases, a maximum size of 4,436,233 bases, a mean size of 142,646 bases, and N50 of 623,820 bases, an N90 of 140,742 and 3.54% N.

Assembled scaffolds were joined using 100 base read data from the sequencing of 11,520 BAC clones. Unique BAC clone reads were extracted and mapped to the assembly scaffolds using megablast (-W 70) [52]. A custom perl script was used to filter out reads mapping at less than 100% of their length and then to merge scaffolds determined to be co-linear. The genetic map of Scaglione and colleagues [15] was used to guide assignment of the vast majority of scaffolds to linkage groups. A few remaining unassigned contigs were assigned to linkage groups using genetic markers from two other sources (Fraser and colleagues [28] and Additional file 12). To enable its use with Red5, markers were first converted to fasta sequences. For Single Nucleotide Polymorphisms (SNP) markers a sequence consisting of the SNP plus 500 upstream and downstream flanking bases was extracted from scaffolds of ‘Hongyang’ [14], grandparent assembly Red Female 1(R Crowhurst, unpublished), an unrelated yellow fleshed A. chinensis assembly CK15_02 (R Crowhurst, unpublished) or the Red5 assembly herein, as appropriate. For each SNP the extracted fasta was named so as to encode the scaffold of origin, the location of the SNP in the scaffold of origin, the sequence region extracted from the scaffold of origin, the linkage group and centimorgan position from the genetic map and the map of origin. The names of the fasta sequences as described are provided in Additional file 12. For EST-based markers [25] the sequence of the EST was obtained from NCBI GenBank. The FASTA sequences for markers were aligned to assembly scaffolds for Red5 using megablast (-W 50) and filtered to remove alignments of less than 98% of overall length and identity before being used for assignment of scaffolds to linkage groups. An ‘all by all’ megablast comparison of scaffolds of ‘Hongyang’ [14] and Red5 was used to assign further Red5 scaffolds to linkage groups using BLAST walking from already assigned Red5 scaffolds. Red5 scaffolds identified as chimeric based on marker evidence (assigned to more than one linkage group or location within a linkage group) were manually inspected, split at identified break points and the component parts re-assigned as supported by evidence (Additional file 12).

To assess the level of DNA sequencing incorporation into the final assembly, the DNA sequencing libraries used as inputs to the anytag software were aligned to the chromosome assemblies using bowtie2 [53] using command line options: --end-to-end --very-fast -I 50 -X 500 --fr --threads 8. To enable comparison with ‘Hongyang’, mapping was repeated using the ‘Hongyang’ chromosomes [14] as the reference.

To assess the accuracy of the assembly, 22 clones from the BAC library of A. chinensis Red Female 1 were selected such that each contained sequence spanning two markers located on Chr25. DNA for each clone was individually barcoded and sequenced using the Life Sciences (Roche) 454 GS-FLX platform. The sequences were assembled using Newbler (version 2.9) and the assembled BAC contigs compared with the Red5 genome assembly using megablast with a word size set at 50. Results were filtered to remove regions of repetitive sequence alignment or alignments with less than 98% match to the Red5 genome sequence. Alignments were converted to GFF3 format and visualised using Geneious (versions 8.1.2) [54]. Additionally, reads from the 9Kb LIMP library were aligned to individual chromosomes for the Red5 whole genome assembly using bowtie2 [53] using the following command line options: --end-to-end --sensitive -k 5 -p 8 --rf -I 1 -X 100000. For each chromosome the distance between mate pairs was visualised using hagfish_blockplot from the software ‘hagfish’ (https://github.com/mfiers/hagfish/). Individual chromosome plots were then cut and pasted to form a montage. Each plot represents the alignment across the entire length of a chromosome. Plots were produced to a standard pixel width irrespective of chromosome length. Green regions indicate mate pairs aligning to the chromosomes within the expected distance for the library. Black indicates regions without mate pair alignment. Pinkish-red indicates regions where the distance between mated paired end reads is shorter (assembly compression relative to physical genome) or greater (assembly expansion relative to physical genome).

Mapping of Red5 RNA-Seq libraries to assembly

The RNA-Seq reads were mapped to the chromosome assemblies of Red5 using the STAR RNA-Seq aligner [55] (version STAR-STAR_2.4.2a) using the command line parameters “--chimSegmentMin 30 --runMode alignReads --alignIntronMin 21 --alignIntronMax 25000 --alignMatesGapMax 25000 --alignEndsType Local”. All RNA-Seq reads were trimmed by 13 bases at the 5´ end prior to use. Reads were additionally trimmed at their 3´ ends when quality score assessment with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) gave quality values below 20.

Gene annotation

Manual curation of gene models was performed as follows. The published gene models for ‘Hongyang’ [14] and A. chinensis ESTs within NCBI GenBank were downloaded and mapped to the assembled pseudo-chromosome sequences of Red5 using GMAP (version 2014–10-22) [31]. The sequences of the Red5 pseudo-chromosomes were repeat masked using RepeatMasker (version open-4.0.5) with options “-e ncbi -pa 30 -s -nolow -species viridiplantae -a -x -poly –gff” and with the RepBase (http://www.girinst.org/) RepeatMasker libraries (20140131). Ab initio gene model prediction was performed using Augustus-3.1 [56] employing command line options “--species = arabidopsis” and evidential hints from both RNA-Seq derived from Red5 (Table 3) and from 47,384 A chinensis EST sequences (downloaded from NCBI GenBank) using described protocols (http://bioinf.uni-greifswald.de/bioinf/wiki/pmwiki.php?n=Augustus.Augustus). The genome sequence, ab initio predicted models, mapped ESTs and ‘Hongyang’ gene models as well as RNA-Seq alignments were imported into WebApollo1 [21] and made available for community-based manual curation via Amazon Cloud. Where evidence suggested multiple isoforms that could clearly be defined each isoform was included. Where it was not possible to differentiate individual isoforms unambiguously a single model was submitted by community annotators. After initial annotation, all models were ported to WebApollo2 (version 2.0.2) (https://github.com/GMOD/Apollo/releases/tag/2.0.2) and each annotation was reviewed manually. This review was followed by computational analysis to identify anomalous annotations such as small (< 15 bp) introns and untranslated exons numbering more than 3. These models were rechecked and often further modified or removed. The manually annotated gene models were named using the prefix Acc (for Actinidia chinensis var. chinensis) and sequentially numbered. The cDNA, CDS, peptide and GFF3 records for each model were committed into a GitHub repository to enable tracking of changes over time. A schema for manual annotation can be found in Additional file 10.

Estimating genome completeness using BUSCO analysis

Benchmarking Universal Single Copy Orthologs (BUSCO analysis) [30] was used to examine genome completeness. Version 2.0 of BUSCO was used with NCBI Blast+ version 2.2.30, AUGUSTUS version 3.2.2 [56] with “--species arabidopsis”, HMMER version 3.1b2 (http://hmmer.org) and the Embryophyta_odb9 dataset (http://busco.ezlab.org/datasets/embryophyta_odb9.tar.gz). BUSCO analysis was also undertaken on the ‘Hongyang’ genome [14].

Comparisons of gene sets

To further evaluate our manually annotated gene models we compared their coding regions with the sequence of 859 clones from cDNA libraries for A. chinensis var. chinensis ‘Hort16A’, within our in-house sequence database which had been previously individually cloned and DNA extracted [25]. The cDNA clones were bi-directionally sequenced in full, using ‘Sanger’ sequencing, and then sequences with greater than 98% similarity were removed using cd-hit-est [57], yielding a comparison set of 812 cDNA sequence with minimum, maximum, and mean lengths of 247, 4506 and 1495 bases respectively and an N50 of 1653 bases (Additional file 7). Coding regions of our manual annotation data set as well as those for ‘Hongyang’ were compared with these cDNA sequences using BLAT (version 36) and GMAP as follows. The longest ORF for each of the 812 cDNAs was extracted and sequences not likely encoding a full length protein based on comparison to NCBI RefSeq Plant (version 76) were removed leaving 550 proteins. Each of the 550 protein sequences was compared to the predicted protein sequences of the 3 gene sets (Red5, original and revised ‘Hongyang’) using BLAT and the best alignments summarised using a custom perl script.

To compare whole gene sets, for each gene set pair (Red5 with original ‘Hongyang’ and Red5 with revised ‘Hongyang’) the following analyses were undertaken: (1) predicted proteins of each pair of gene sets were compared to each other using BLAT, (2) CDS sequences for each gene set pair were compared using BLAT, (3) the CDS were aligned to the genome sequence using GMAP (version 2017–06-20). A custom perl script was then used to summarise these analyses by first seeking the best protein:protein alignments, then the best CDS:CDS alignments and finally the alignment of the query to the whole genome sequence. The summarisation perl script takes into account query and target sequence lengths for protein:protein and CDS:CDS alignments as well as genome alignment co-ordinates from the GMAP alignments in order to yield the metrics for best alignments including: number of alignments assigned to bins based on percentage identity, the number showing a length variance (query equal/shorter/longer than target), number of queries encapsulated with the UTR of a target gene model, number of queries overlapping the 5 or 3 UTR of a target gene model and the number missing from results files for these analyses. Additional file 13 summarises the raw results for the best BLAT gene set match to each cDNA protein. BLAT results presented in Additional file 13 were extracted using “awk ‘{Percent=($1/$11)*100; print $10''\t''$11'' \t''$1''\t'' Percent''\t''$14''\t''$15}’” then the difference in length between summary data columns 2 and 6 added in Microsoft Excel.

Phylogenetic analysis

Gene models were selected and aligned in Geneious (R10.0.3) (www.geneious.com) using Geneious Alignment (with free end gaps), Gap opening penalty 30, extension penalty 0 and refinement iterations 2. Alignable regions were extracted, realigned and clustered using PHYML [58] default settings. Data from 1000 bootstrap sets are presented.