A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds

Contig generation with PacBio and 10x chromium

We generated data with PacBio, 10x Chromium, BioNano, and Hi-C. The PacBio and 10x Chromium sequences were first used to produce separate independent assemblies using FALCON and Supernova respectively (see Methods). The PacBio assembly had the highest contiguity of these single-technology assemblies, with 429 primary contigs of average size 0.520 Mbp and N50 of 3.09 Mbp (Table 1). We next scaffolded the PacBio assembly with 10x data using ARCS [35] and LINKS [36], and oriented contigs and scaffolds on a genetic map followed by additional gap filling with PBJelly [37]. The contiguity of this assembly version (Amel_HAv1) was significantly improved compared to both the individual 10x and PacBio assemblies. The longest Amel_HAv1 contig is 13.4 Mbp, 40 times longer than in the longest contig in Amel_4.5. N50 of the HAv1 is 5.167 Mbp, compared to 0.046 Mbp for Amel_4.5 (Fig. 1; Table 1).

		Amel_4.5	10 × a	PacBio	Amel_HAv1c	Amel_HAv2d	Amel_HAv3e
Size	Total (Mbp)	229.12	217.80	223.24	225.21	225.23	223.86
Contigs (all)	N	16,501	20,240	429	330	331	228
	Longest (Mbp)	0.333	0.288	9.726	13.399	13.399	13.400
	Mean (Mbp)	0.014	0.011	0.520	0.682	0.684	0.974
	N50 (Mbp)	0.046	0.031	3.086	5.167	5.167	5.381
	L50 (n)	1390	1968	23	14	14	13
Scaffolds (all)b	N	5644	9734	–	–	280	177
	Longest (Mbp)	4.736	3.297	–	–	27.79	27.77
	Mean (Mbp)	0.041	0.024	–	–	0.816	3.340
	N50 (Mbp)	0.997	0.589	–	–	11.33	13.62
	L50 (n)	65	116	–	–	8	7
Scaffolds (anchored to nuclear chr.)	N	340	–	–	–	36	16
	Longest (Mbp)	4.736	–	–	–	27.79	27.77
	Mean (Mbp)	0.598	–	–	–	6.21	13.79
	N50 (Mbp)	1.209	–	–	–	11.60	13.62
	L50 (n)	52	–	–	–	7	7

^aLinked-read sequences taken as scaffolds. Contigs derived from splitting scaffolds on Ns

^bIndividual unplaced fragments counted as scaffolds

^cPacBio+10x

^dPacBio+10x + BioNano

^ePacBio+10x + BioNano+Hi-C

Scaffolding with BioNano and hi-C

We performed scaffolding of the Amel_HAv1 contigs using BioNano data to produce version Amel_HAv2. This version contains 26 hybrid scaffolds with N50 of 11.3 Mbp and the longest scaffold of 27.8 Mbp. In total 96 out of 171 available BioNano genomic maps could be used to scaffold contigs. The remainder could not be anchored to contigs, or did not link multiple contigs. A total of 77 out of 328 Amel_HAv1 contigs were scaffolded using the BioNano genomic maps, whereas the remaining contigs were retained unchanged. Six of the sixteen chromosomes were recovered as single scaffolds and each chromosome was represented by an average of 2.2 scaffolds. The remaining unplaced contigs were comparatively short and highly repetitive.

We conducted additional scaffolding using the genetic map AmelMap3 [38] and Hi-C data, followed by gap filling and polishing in order to produce version Amel_HAv3. In this final version, each chromosome is represented by a single scaffold, comprised of an average of 4.2 contigs. Chromosomes 4 (13.4 Mbp) and 15 (9.5 Mbp) are recovered as single contigs, including the distal telomeres (see below). For comparison, in Amel_4.5, the chromosomes are comprised of 340 anchored scaffolds. Contigs are named after linkage group and order on the genetic map, i.e. Group1_2 for the second contig on linkage group 1. A full list of scaffolds, contigs and their length and placements is provided in Additional file 1: Table S1. A visual overview of the 16 chromosomes is presented in Fig. 2.

Congruence of assembly with the genetic map

The order of genetic map markers in the linkage map AmelMap3 [38] was compared to their order on the Amel_HAv3 chromosome-length scaffolds. Out of a set of 4016 paired primer sequences for 2008 microsatellite markers (Additional file 1: Table S2), we found that 301 primers for 268 markers did not map to the assembly (7.5% of primers; 13.3% of markers), including both primers for 33 markers. Thus 1975 marker loci (98.4%) could be positioned along the chromosomes (avg. 123 markers per chromosome). Out of these, 1885 (95.4%) are congruent and collinear between Amel_HAv3 and the genetic map and the scaffolds are nearly fully consistent with the order of contigs suggested by the genetic map (Additional file 1: Table S3). We find a small fraction (0.9%) of the markers to be ambiguous. The primer pairs were originally designed to amplify polymorphic microsatellites and are expected to map close together on the chromosomes and not overlap with other pairs. The BLAST targets were > 1 kbp apart for only 10 primer pairs (0.5%) and for 8 pairs (0.4%) they were overlapping.

However, we also detected minor unresolved incongruences inside or between adjacent contigs. A total of 72 markers (3.6%) have inconsistent placements in Amel_HAv3. These include cases where a small number of adjacent markers were locally arranged in the opposite physical order along contigs, compared to the expected order in the genetic map or where markers from different adjacent contigs were mixed at their borders, producing interleaved or nested contigs with respect to their order in the genetic map. Removing markers at zero genetic distances to their adjacent markers (n = 241) reduced this rate of inconsistency to 2.5%, suggesting that the original order of some of these markers in the genetic map is itself ambiguous. Interleaved/nested contigs were observed within 5 chromosomes: the 0.4Mbp contig Group6_2 appears to be partly discontinuous and nested within Group6_1 on chromosome 6; contig Group7_2 overlaps with the end of Group7_1 on chromosome 7; a single-marker from Group10_6 is associated with Group10_5 on chromosome 10; Group12_1 and Group12_2 are interleaved across a 0.1–0.2 Mbp region on chromosome 12; and a ~ 0.3Mbp segment of Group13_5 is found within Group13_6. These inconsistencies and marker primers that could not be placed on the new assembly may indicate unresolved assembly errors or other sequence differences around these microsatellite loci (e.g. missing or divergent target sequence between this assembly and that used to produce the markers). Alternatively, they may reflect natural structural variation between the sample used for this assembly and those used to produce the genetic map.

Comparisons of anchored and unplaced contigs in Amel_4.5 and Amel_HAv3

In Amel_4.5, 16.7 Mbp (7.3%) of the sequence is marked as repetitive and unplaced contigs have higher levels of repeat sequence than chromosome-anchored contigs (Table 2). In Amel_HAv3, the overall amount and proportion of repeats has increased to 17.4 Mbp and 7.9%. In comparison to the overall addition of sequence to chromosomes (+ 10%; 219.4 in Amel_HAv3 vs. 199.7 Mbp in Amel_4.5), we find that repeat sequence has been added at twice this proportion (+ 21%; 16.5 Mbp vs. 13.6 Mbp), indicating that we have incorporated sequence with higher levels of repeats than the genomic background into chromosomes.

Several features distinguish contigs that we were unable to incorporate into the genetic map or scaffolds (Table 2). These contigs are lower in GC content, have a larger proportion of repetitive sequence and have lower mappability. These features are also present in Amel_4.5, but are more pronounced in Amel_HAv3. For instance, repeat content is 2.9-fold higher among the unplaced vs. anchored Amel_HAv3 contigs compared to 1.54-fold higher in Amel_4.5 (Table 2). These repeat sequences remain difficult to place even with current long-read technologies.

BUSCO gene content

We compared the respective completeness of the Amel_4.5 and Amel_HAv3 assemblies by counting the number of universal single-copy orthologues detected in either assembly with BUSCO [39]. Overall, Amel_HAv3 has a slightly larger number of BUSCO genes compared to Amel_4.5 (Table 2). However, in Amel_4.5, 5–6% of these BUSCOs are detected among unplaced contigs, whereas only 0.4–0.8% of these occur in unplaced contigs in Amel_HAv3 (Additional file 1: Table S4). The hybrid assembly therefore represents a significant improvement in terms of the proportion of conserved genes located in genome scaffolds.

The mitochondrial genome

We recovered a complete mitochondrial genome (16,463 bp) and could detect and label all features along the sequence (13 coding genes; 22 tRNAs; 2 rRNAs) using a combination of BLAST [40] and MITOS [41]. All coding genes and rRNAs, and most tRNAs (n = 15), were accurately detected using BLAST (< 6 bp missing from canonical models). All tRNAs were detected near-full length using MITOS (< 3 bp missing from canonical models). All features were found to be in full synteny with previous assemblies [34, 42]. The Amel_HAv3 mitochondrial sequence is 120 bp longer than in these assemblies. After aligning the sequences, we found that most of the length difference is explained by three major intergenic indels: i) a 16 bp deletion between COX3 and tRNA-Gly; ii) a 190 bp hyper-repetitive insertion in the AT-rich region (%AT = 96.9) next to the small ribosomal subunit; and iii) a 39 bp deletion in the same region. The remaining 15 bp are due to small scattered 1–3 bp indels. The 190 bp insertion was likely not possible to assemble before with Sanger or short-reads. The mitochondrial genome and structural variants are presented in Additional file 2: Figure S2 and feature coordinates are provided in Additional file 1: Table S5.

Repeat content

The honeybee genome has relatively few repeats compared to other insects (8%; Table 2). In both this and the previous assemblies (Amel_HAv3 and Amel_4.5), we find 12.8 Mbp of simple repeats/low complexity regions with RepeatMasker, representing 5.6% of the overall sequence and about 75% of all repeat-masked output (Additional file 1: Table S6). The remaining share (5 Mbp) consists of longer interspersed DNA transposons, long/short interspersed nuclear elements (LINE/SINEs), long terminal repeats (LTRs), RNA sequences and other minor repeat classes. In agreement with previous analyses of transposable elements in honeybees [34], we find that DNA transposons are the major repeat class (3.1 Mbp; 66% of all interspersed repeats; 1.4% of the assembly; Fig. 3a; Additional file 1: Table S6), and that mariner transposons are the most common element within this class (1.74 Mbp; 56% of DNA transposons). Many repeats occur at approximately the same frequency in both assemblies under our analytical conditions (Fig. 3b-c), although some repeat classes occupy larger proportions of the genome. For instance, DNA transposons are only 1.02 times more frequent but occupy 1.25 times more space in Amel_HAv3 compared to Amel_4.5. Likewise, rRNA sequences occupy over two times as much sequence but occur at nearly the same frequency (Additional file 1: Table S6). This discrepancy suggests that many repeat motifs are individually longer in Amel_HAv3 than in Amel_4.5.

The most striking difference in repeat annotation in Amel_HAv3 is the addition of a large number of AvaI (547 bp; n = 229) and AluI (176 bp; n = 1315) repeats (Fig. 3b; Additional file 1: Table S7). These repeats have previously been estimated to represent 1–2% the honeybee genome using Southern blotting and FISH, and to be clustered close to centromeres (AvaI) and the short-arm telomeres (AluI) [30, 43]. We detect 6.5 times more AluI repeat sequence in Amel_HAv3 than in Amel_4.5 and 11 times more AvaI sequence (Fig. 3c), although we are unable to fully assemble and map the complete sets because many of the repeats occur in unplaced contigs (89% of AluI and 41% of AvaI repeats, respectively). The enrichment is lower by fragment count rather than overall sequence length (5.2-fold for AluI and 5.1-fold for AvaI). This is likely explained by higher repeat fragmentation in Amel_4.5, inflating repeat counts: only 30% of AluI repeat matches are > 160 bp in Amel_4.5, compared to 78% in Amel_HAv3 (Additional file 3: Figure S3A). Likewise, average divergence from the canonical AluI repeat is 15% in Amel_4.5 but only 3.9% in Amel_HAv3 (Additional file 3: Figure S3B). For the AvaI repeats, only 2% are > 500 bp in Amel_4.5 vs. 73% in Amel_HAv3, and divergence is 21% vs. 6.6% (Additional file 3: Figure S3C-D).

In Amel_HAv3, we find that AluI and AvaI repeats tend to cluster into extended tandem arrays (see Fig. 2 for their distribution on anchored contigs), often without any extra bases inserted between copies. The longest such anchored gap-free AluI array occurs on contig Group2_1 at the start of chromosome 2 and spans 80 adjacent full-length copies reiterated across 14.0 kbp (89 AluI repeats occur in the region; Fig. 4a). Half of all AluI repeats occur in tandem arrays of at least 21 copies (~ 3.7 kbp). Two of these are at the short-arm ends of scaffolds (chrs. 2 and 12; Fig. 2), while the rest are on unplaced contigs, indicating that these extensive AluI repeats are associated with the short-arm telomeres, as suggested previously [30, 43]. In comparison, the longest contiguous AluI region in Amel_4.5 spans only 9 repeats (~ 1.4 kbp) on an unplaced contig. Likewise, the longest AvaI region in Amel_HAv3 occurs on chromosome 1 (contig Group1_3) and spans 26 copies across 14.1 kbp (Fig. 4b), and 50% of AvaI repeats occur in tandem arrays of at least 9 repeats (~ 4.4 kbp). No such gap-free AvaI array is detected in Amel_4.5. These improvements in Amel_HAv3 underscore the major advantage that long reads have over previous short-read technologies in resolving and representing highly repeated sequences.

Alignment to previous builds

To further characterize differences between the genome assemblies, we aligned Amel_4.5 sequences against Amel_HAv3 using Satsuma [44]. Overall, alignments were produced against 94.3% of Amel_HAv3 (see Additional file 3:. Figure S3 for full alignment maps between assemblies). Chromosomal sequence was more frequently aligned (95.4% of 219.4 Mbp in Amel_HAv3) than unplaced contigs (44.4% of 4.45 Mbp in Amel_HAv3), which is consistent with these relatively repetitive contigs containing sequence that is not well represented in Amel_4.5. For sequences that had been associated with chromosomes in both assemblies, we found that 191.6 Mbp of alignments originated from the same chromosome in either assembly (99.4% of chromosome-to-chromosome alignments; 86% of all Amel_HAv3), while only a small fraction (1.23 Mbp; 0.56% of Amel_HAv3) originated from different chromosomes (Fig. 5a), suggesting largely consistent mapping of data. About 16.4 Mbp of sequence that had previously been unplaced in Amel_4.5 now aligned against Amel_HAv3 chromosomes, corresponding to 7.5% of the total Amel_HAv3 assembly (Fig. 5a). For comparison, we found that the opposite pattern was very uncommon: only 0.148 Mbp of alignments was mapped to chromosomes in Amel_4.5 but is unplaced in Amel_HAv3 (0.07% of Amel_HAv3). About 10.3 Mbp (4.7% of Amel_HAv3) was anchored to chromosomes but had no matching sequence in Amel_4.5 (conversely, 6.6Mbp of chromosomal sequence in Amel_4.5 is not matched in Amel_HAv3). Alignments were produced for 1.98 Mbp of contigs that are unplaced in both assemblies (0.9% of Amel_HAv3), whereas 2.32 Mbp unplaced Amel_HAv3 contigs did not align against Amel_4.5 (1.1% of Amel_HAv3; Fig. 5a).

Aligned sequences that are anchored to the same chromosomes in either assembly have the highest average mappability scores (0.994) and GC content (34.1%; Fig. 5b), characteristic for high-complexity/low-repeat sequence that is most amenable to assembly via last-generation technologies. Sequence that has been incorporated into chromosomes only in Amel_HAv3, but is unplaced in Amel_4.5 or unmatched/unaligned with Amel_4.5 sequence, has significantly lower GC-content, and in the latter case also lower mappability. Both aligned and unaligned sequences that we are unable to place on chromosomes have reduced mappabilities and GC-content compared to genome genomic background (Fig. 5b). Sequence that has switched chromosomes between assemblies has intermediate values for these statistics.

We find that newly anchored sequences or sequences that are still unplaced are significantly enriched for both simple and interspersed repeats (Fig. 5b; see Additional file 4: Figure S5 for the density of individual repeat classes). Chromosomal regions built from sequences that were unplaced in Amel_4.5 or unaligned to Amel_4.5 sequence represents 12% of the genome but contain 17% of simple repeats, 42% of DNA transposons, 25% of LTRs, 35% of satellites and 59% of AvaI repeats (Fig. 5a). Regional occurrence and enrichments for repeat-classes and their sub-classes can be found in Additional file 1: Table S8.

Distal telomeres

The telomeric repeat motif TTAGG is expected to occur as tandem arrays at the tip of the distal long-arm telomeres of all honeybee chromosomes. Distal telomeres were previously characterized from relatively short and fragmented sequences spanning only a few hundred base pairs of TTAGG repeats at the tips of five long-arm chromosomes in assembly Amel_4.0 [33], but manual scaffolding connected them to all but one long-arm chromosome tip [45]. We scanned for TTAGGs across 10 kbp windows in Amel_HAv3 and detected large clusters (on average 1177 repeats; ~ 5.7 kbp) at the very ends of the long arms of 14 chromosomes (all except chromosomes 5 and 11; Fig. 2; Additional file 1: Table S1). While TTAGG/CCTAAs are rare across the genome (about 8 motifs per 10 kbp or ~ 0.4% of the genomic background; Fig. 2), the outermost 1–2 windows of these chromosomes contain on average 1043 motifs per 10 kbp (52% of the sequence; 130-fold enrichment; Fig. 2). The longest telomeric repeat region was assembled for chromosomes 3 and 8, containing 2142 and 1994 copies of the motif, respectively. For the metacentric chromosome 1, we detected TTAGG repeats at both ends of the chromosome (the reverse complement motif CCTAA at the start of the chromosome), which is consistent with the hypothesis that this large chromosome has formed from fusion of two acrocentric chromosomes and harbors two distal telomeres [33, 45].

We extracted and aligned the sequences of all distal telomeres with TTAGG arrays using MAFFT (n = 15, including both telomeres on chromosome 1), including ~ 4 kbp of the upstream subtelomeric region, and scanned the sequences for shared properties. Taking the sequence at chromosome 8 as reference, we find that the first 2kbp downstream of the start of the telomere is enriched for TCAGG, CTGGG and TTGGG variants (Fig. 6a, c). These polymorphisms are gradually replaced by the canonical TTAGG repeat moving towards the distal ends of the telomeres, where the average pairwise divergence between telomeres accordingly is much reduced: from 12% at < 2 kbp away from telomere start to 2.4% at 2–4 kbp away (Fig. 6d). We recover a relatively conserved 3 kbp subtelomeric region upstream of the junction (avg. pairwise divergence 14%; Fig. 6d). The subtelomeres contain two larger shared motifs just upstream of the junction telomere junction (Fig. 6): i) a ~ 350 bp (213-520 bp) fragment is located 100 bp upstream of the junction and has moderate similarities towards a 4.5kbp LINE/CR1 retrotransposon originally characterized in Helobdella robusta (CR1-18_HRo; avg. ~ 74% identity; ii) a highly conserved and GC-rich 400 bp sequence (avg. div. 5.5%) is located further upstream but does not have significant similarities with any sequences in RepeatMasker or NCBI GenBank. Chromosomes 5 and 11 do not contain arrays of TTAGGs in Amel_HAv3, but terminate with subtelomere sequences that include the conserved motif (identified with BLAST). Three unplaced contigs contain a large number of motifs: the 18 kbp GroupUN_199 has 1177 TTAGGs, the 16 kbp GroupUN_7 has 909 CCTAAs and GroupUN_198 has 82 CCTAAs. (Additional file 1: Table S1). No other 10 kbp window contains > 30 such motifs among the unplaced contigs. It is possible that these three contigs belong to the truncated chromosomes. Both GroupUN_7 and GroupUN_198 associate with chromosome 11 in the Hi-C dataset. GroupUN_198 also contains a > 2.6 kbp subtelomeric subsequence (labeled with BLAST). Mate-pair reads with TTAGGs have previously been linked to the tip of chromosome 11 [45].

For comparison, we scanned Amel_4.5 for TTAGGs and subtelomeric sequence to locate telomeres in this assembly. In Amel_4.5, we find subtelomeres on short contigs (average length of 27kbp) located at the tips of the outermost scaffolds of 13 chromosomes. We detected TTAGG clusters with on average 34 motifs per telomere near five of these, whereas the rest had repeat densities that were indistinguishable from background levels. The distal telomere sequences in Amel_HAv3 are therefore 35 times longer (1177/34) than those in Amel_4.5.

Centromeres and proximal telomeres

AvaI and AluI repeats have previously been suggested to indicate the positions of centromeres and proximal short-arm telomeres, respectively, in the honeybee genome [30, 43]. Although many AvaI and AluI repeats remain unmapped (see above), we find that the mapped repeats cluster toward the tips of the short-arms of most acrocentric (2–6, 9, 14–15) and the center of metacentric chromosome 1 and possibly submetacentric chromosome 11 (Fig. 2; Additional file 1: Table S7). These locations are largely similar to previous FISH labeling of these sequences for all but two chromosomes (10 and 16; Fig. 2). We find that short-arm telomeric AluI repeats often co-occur with centromeric AvaI repeats (e.g. chromosomes 2, 4, 5 and 6). This is consistent with fluorescent labeling that also suggests that the proximal telomeres blend with centromeres in many acrocentric chromosomes [30]. However, the assembly often terminates at or near these clusters, sometimes before reaching into the proximal telomere (e.g. chromosomes 14 and 15; Fig. 2).

The distribution of the putatively centromeric AvaI repeats in Amel_HAv3 overlaps or co-occurs with experimental mapping of centromeres from patterns of recombination and heterozygosity in half-tetrads of the clonal Cape honeybee A. m. capensis (e.g. chromosomes 2, 4, 11; Fig. 2) [31]. The high contiguity in Amel_HAv3 now facilitates further characterization of the putative centromeric regions. All mapped AvaI clusters with more than two repeats (n = 11; Fig. 2) are embedded in megabase-scale regions with reduced GC content compared to the rest of the genome (22.7% vs. 34.6%; average length 2.3 Mbp; delineated by 100kbp windows with GC < 32.7%; Fig. 2). Sequences up to 1–2 Mbp away from the AvaI clusters have significantly reduced %GC and increased density of simple repeats and DNA transposons, compared to the genomic background (> 2 Mbp away; p < 0.05; 2000 bootstrap replicates of data intervals; Fig. 7). Patterns of centromeric enrichment were unclear for the rarer repeat classes. Similar low-GC blocks were detected in chromosomes 13 and 16, although only a single or no AvaI repeat, respectively, was mapped to these regions. The low-GC centromere-associated regions together span 42 Mbp of the genome and are among those that appear to have been particularly poorly assembled before: these regions constitute 19.3% of the genome but contain 38% of all sequence that is unmatched against Amel_4.5 and 95% of all sequence that was unplaced in Amel_4.5 (Additional file 5: Figure S4A). These regions have more than doubled in size compared to Amel_4.5.

We next used the genetic distances previously inferred between the genetic map markers to compare recombination rates inside and outside of these regions. Across the genome, we estimate the average recombination rate to be 21.6 cM/Mbp (n = 1735 congruent marker pairs), close to what has been estimated before in honeybee [38, 46, 47]. Compared to these background levels, recombination rates are significantly reduced across both sets of centromere mappings: to 14.6 cM/Mbp in the half-tetrad experiment from [31] and to 7.9 cM/Mbp from our assessment of AvaI and GC-content (Fig. 8). In contrast to the FISH results, we also detect several small AluI clusters close to long-arm telomeres (Fig. 2). However, compared to the repeats at the proximal telomeres, these hits are fewer (32 vs. 130), shorter (106 bp vs. 162 bp) and more divergent (16% vs. 4.5%) on average (Additional file 3: Figure S3A-B), which could indicate excess spurious hits or degenerate elements.

We do not find TTAGGs associated with proximal telomeres, suggesting they are either not present at the short-arms of the honeybee chromosomes or only occur in unmappable sequence. To address this, we manually inspected mate-pairs sequenced from decade-old fosmid libraries that were prepared for the original assembly [33]. Fosmid reads containing AluI repeats were found to likely have AluI mate pairs, indicating very long strings of AluIs that supersede the length of the arrays in the hybrid assembly. Interestingly, out of 19 mate pairs containing the TTAGG motifs and linked back to telomere regions, only 7 anchored to distal telomeres, while 12 contained AluI repeats. This suggests they belong to proximal telomeres, although no individual read contained directly observable junctions between AluI and TTAGG motifs. This independent evidence nevertheless suggests the presence of TTAGG repeats beyond the currently mappable regions of AluI repeats on the short-arm telomeres. Moreover, the CCTAA enriched unplaced contig GroupUN_7 (see above) also contains 28 AluI repeats, and could potentially be a proximal or mis-joined contig. Because our assembly of these regions between the centromeres and the short-arm telomeres remains incomplete, most of the unplaced contigs are inferred to belong in these regions.

Library preparation and data production

We produced data using Pacific Biosciences SMRT sequencing (PacBio), 10x Chromium linked-read sequencing (10x), BioNano Genomics Irys optical mapping (BioNano) and a Hi-C chromatin interaction map (Phase Genomics). DNA extracted from a single drone pupa from the DH4 line was used for the first three of these methods (a different drone for each method). These individuals were brothers of the individuals from the DH4 line used for previous honeybee genome assembly builds [33, 34]. The sample used for Hi-C was an individual from an unrelated managed colony with a similar genetic background as the DH4 line (mixed European) collected from the USDA-ARS Bee Research Laboratory research apiary.

To prepare DNA for the PacBio and 10x sequencing, we first lysed cells from 20 to 120 mg of insect tissue. This was done by grinding in liquid nitrogen followed by incubation at 55 °C in cell lysis solution (25 ml 1 M Tris-HCl, pH 8.0; 50 ml 0.5 M EDTA, pH 8.0; 0.5 ml 5 M NaCl; 12.5 ml 10% SDS; 162 ml molecular grade water) and proteinase K. The solution was then treated with RNase A. Proteins were then precipitated using Protein Precipitation Solution (Qiagen) and centrifugation at 4 °C. DNA was precipitated from the resulting supernatant by adding isopropanol and ethanol and centrifugation at 4 °C.

We generated a 10 kb PacBio library that was size-selected with 7.5 kb cut-off following the standard SMRT bell construction protocol according to manufacturers recommended protocols. The library was sequenced on 29 SMRT cells of the RSII instrument using the P6-C4 chemistry, which generated 10.2 Gb of filtered data. N50 subread length was 8.8 kb. A 10x GEM library was constructed from high-molecular-weight DNA according to manufacturers recommended protocols. The resulting library was quantitated by qPCR and sequenced on one lane of a HiSeq 2500 using a HiSeq Rapid SBS sequencing kit version 2 to produce 150 bp paired-end sequences. This resulted in 127,440,953 read pairs (38Gb of raw data).

High-molecular-weight DNA was extracted in situ in agarose plugs from a single drone pupa following BioNano Genomics guidelines. Plugs were cast and processed according to the IrysPrep Reagent Kit protocol with the following specifications and modifications; a 7-day proteinase K treatment in lysis buffer adjusted to pH 9.0 with 2 μl BME per ml buffer. The BspQI NLRS reaction was processed according to protocol, stained overnight and immediately loaded on 2 flow cells for separation on the BioNano Irys system. In total 1,214,651 molecules were scanned with N50 of 210 kbp.

DNA for the Hi-C experiment was prepared at Phase Genomics. The sample was incubated at 27 °C for 30 min with periodic mixing by inversion. Glycine was added (final concentration of 0.1 g/10 mL) to quench crosslinking. After an additional incubate at 27 °C for 20 min with periodic inversion, the sample was pelleted by centrifugation, the supernatant was removed and the sample was kept at − 20 °C prior to processing. This procedure results in extraction of native cross-linked chromosomes from bee cells, disruption using endonucleases and linking of adjacent strands via biotinylated junctions. The samples were then sequenced on an Illumina HiSeq instrument.

Assembly pipeline

In order to determine the best way to utilize the data, we first generated assemblies using the PacBio and 10x Chromium data independently (see below). As the PacBio assembly had far superior contiguity, we designed a pipeline to begin with this assembly and then use the 10x linked-reads to connect and combine contigs. These contigs were then scaffolded using the BioNano optical map data, with additional checks for consistency with Hi-C data and a genetic map [38]. Gap filling and polishing steps were also included to maximize contiguity and accuracy. Full details of the pipeline are presented below and summarized in Additional file 6: Figure S1.

We imported PacBio raw data into the SMRT Analysis software suite (v2.3.0) (Pacific Biosciences, CA) and generated subreads. All sequences shorter than 500 bp or with a quality (QV) < 80 were filtered out. The resulting set of subreads was then used for de novo assembly with FALCON v0.5.0 [65] using pre-assembly length cutoff of 7 kbp. Since the genomic DNA originated from haploid drone we kept only primary contigs generated by FALCON and removed 14 contigs shorter than 2kbp before further analysis. The resulting set of contigs was polished twice using Quiver via SMRT Analysis Resequencing protocol [65]. The resulting PacBio assembly consisted of 429 contigs with N50 of 3.1Mbp and largest contig being 9.7Mbp.

To create the 10x Chromium assembly we used Supernova 1.1.4 on the 10x Chromium linked read data [66] with default parameters. The resulting assembly had 9734 scaffolds with N50 of 0.59Mbp and the longest scaffold was 3.2Mbp. This assembly was not used to create the final assembly and we instead used the 10x linked-read data to extend the PacBio assembly generated in the previous step. We ran the ARCS+LINKS Pipeline [35, 36] to utilize the barcoding information contained in 10x linked reads. First, we mapped 10x reads to PacBio contigs using LongRanger 2.1.2 (10X Genomics, CA). ARCS v1.0.1 was then used to identify pairs of contigs with evidence that they are connected based on the observation of linked reads from the same molecule. Default parameters were used, except for modifying barcode read frequency range (−m 20–10,000). The results of ARCS were processed with the LINKS v1.8.5 scaffolding algorithm to constructs scaffolds based on 10x read pairing information. We adjusted the –a parameter, which controls the ratio of barcode links between two most supported graph edges, to 0.9. The ARCS+LINKS pipeline produced 299 scaffolds with N50 of 8.8Mbp and longest scaffold of 13.3Mbp.

We compared the PacBio+10x assembly to the genetic linkage map AmelMap3 [38] by determining the position of 2008 microsatellite markers using BLAST [40, 67]. This enabled us to assign 49 scaffolds to one of 16 linkage groups. The remaining sequences were designated as unplaced. Furthermore, we used genetic map information to order, orientate and to join adjacent scaffolds belonging to the same linkage group by introducing arbitrary gap of 2000 Ns. We then used PBJelly from PBSuite v15.8.24 [37] to perform a first round of gap filling using all PacBio reads. PBJelly closed 87 (67%) gaps within scaffolds due to joins made by ARCS+LINKS and 16 (48%) of the gaps that were introduced between adjacent scaffolds on the basis of proximity according to the genetic map. In order to minimize possibility of freezing scaffolding errors we then split scaffolds on remaining gaps. This stage of assembly resulted in assembly version Amel_HAv1 which had 330 contigs with N50 of 5.6Mbp and longest contig of 13.4Mbp.

The BioNano raw data were assembled using the BioNano Solve (v3.1.0) assembly pipeline (BioNano Genomics, CA) on a Xeon Phi server resulting in 171 genome maps with cumulative length of 285 Mbp and N50 of 2.2 Mbp. This data set was combined with Amel_HAv1 by running the BioNano Solve v3.1.0 hybrid scaffolding pipeline. The BioNano software identified 7 conflicts between optical maps and Amel_HAv1. All of the conflicts could be traced back to original FALCON assembly and were confirmed to be chimeric. Therefore, we chose to resolve these conflicts in favor of the BioNano optical maps. The resulting hybrid assembly had N50 of 11.3Mbp and a longest scaffold of 27.7Mbp length. This version of the assembly was designated Amel_HAv2.

Hi-C read pairs were aligned to the initial assembly using BWA-MEM [68] with the “-5” option. Unmapped and non-primary alignments were excluded using samtools [69] with the “-F 2316” filter. We next performed scaffolding with Hi-C data using the Proximo pipeline (Phase Genomics), which builds on the LACHESIS scaffolding package [14]. In total 149 out of 280 Amel_HAv2 scaffolds could be grouped into 16 clusters out of which 14 clusters contained scaffolds that were previously assigned to linkage groups and each had cumulative size of predicted chromosome length. Two clusters contained short contigs without linkage group assignments. Two large Amel_HAv2 scaffolds were not a part of any Hi-C cluster, most likely due to fact that they already represented complete chromosomes. Overall, we observed completely accurate assignment of scaffolds to chromosomes based on comparison with genetic map. However, there were a number of errors in orientation and order of scaffolds within a Hi-C cluster. In total there were 22 such errors affecting 8 Hi-C clusters. Therefore, in the final chromosome scale scaffolds we ordered and orientated Amel_HAv2 scaffolds based on genetic map information. Scaffolds from the same linkage group and same Hi-C cluster were joined by introducing an arbitrary gap of 200 Ns in the orientation and order indicated by the genetic map.

We performed an additional round of gap filling using PBJelly to fill gaps generated by the previous scaffolding steps using all PacBio reads, which closed 12 (20%) additional gaps within the scaffolds. New sequences introduced during gap filling steps originate from non-error corrected subreads. To remove potential sequencing errors, the whole assembly was once more subject to two rounds of Quiver polishing. Subsequently 89 unplaced contigs that had fewer than 50% of their bases polished were removed from the final assembly. In addition, two unplaced contigs were identified as mitochondrial and removed. The final assembly consisted of 16 chromosomal scaffolds with a total of 51 gaps, 160 unplaced contigs and a mitochondrial sequence. This final data set was designated as hybrid assembly Amel_HAv3.

Assembly characterization and analysis

After the scaffolding was completed, the congruence between the genetic map and the final assembly was reassessed. The primer sequences for the microsatellite markers in AmelMap3 were again fitted against the genome using BLAST. The physical positions and order of the markers along and between contigs was compared to their expected order in the linkage map.

The assembly was scored for base composition, mappability and repeat content. These metrics may correlate with sequences that are challenging to assemble. We compared the chromosome-anchored and unplaced sequences of the published reference assembly (Amel_4.5; [34]) to the new hybrid assembly (Amel_HAv3) for these properties. We computed average GC content across whole assemblies, arbitrary regions and along non-overlapping windows of different sizes (1 kbp; 10 kbp). We then used GEM v1.315b [70] to compute the mappability (or uniqueness) of short non-degenerate (0 bp mismatch) 50 bp kmers across the assemblies. Every base is annotated for an average mappability score computed from overlapping kmers. We computed average mappability scores across windows as above.

We used BUSCO v2.0.1 [39] to compare the completeness of Amel_4.5 and Amel_HAv3 by assessing the number of expected and detected single-copy orthologs in either assembly as inferred from the OrthoDB v9.1 [71]. Two core sets of BUSCOs (near-universal single-copy orthologs) were used: Metazoa (n = 978 BUSCOs) and Hymenoptera (n = 4415 BUSCOs).

We used RepeatMasker v4.0.7 [72] to annotate simple and interspersed repeat content. We deployed the RMBLAST-NCBI search engine to scan for animal repeats (−species metazoa) in the 20,170,127 release of the Repbase database [73]. The query and database was extended (−lib) to include the consensus motifs of two tandem repeats associated with centromeres (AvaI; 547 bp; X89539) or proximal short-arm telomeres (AluI; 176 bp; X57427), respectively, of honeybee chromosomes [30, 43]. These elements were named after the bacterial restriction endonucleases originally used to detect them (AvaIR from Anabaena variabilis; AluIR from Arthrobacter luteus) but are unrelated to the similarly named Ava and Alu SINE class repeats of other taxa and not detected using the Repbase database. The canonical AvaI repeat consists of four highly similar sub-repeats, resulting in spurious overlapping annotations when AvaI repeats occur in tandem. We therefore parsed the pre-ProcessRepeats output (ori.out-file) separately to extract non-overlapping AvaI repeats. Simple repeats and low-complexity sequence as annotated by RepeatMasker were considered together.

In order to locate distal telomeres, we estimated the density of the short telomeric repeat motif TTAGG/CCTAA [45] across 10kbp windows in both Amel_HAv3 and Amel_4.5 using a custom Perl script. Distal telomeric and subtelomeric regions (<5kbp upstream of putative telomere start) were then extracted and aligned with the L-INS-i algorithm in MAFFT v7.310 [74]. The sequence in chromosome 8, which has among the longest telomeric repeat regions, was taken as a profile and columns with gaps in this sequence was removed from the alignment. The average pairwise sequence divergence, GC-content and density of TTAGG and variant repeats (TCAGG, CTGGG, TTGGG) was then estimated across 100 bp windows along the alignment and 95% confidence intervals were computed by bootstrapping the sequences (n = 2000 replicates). We searched Amel_4.5 for a conserved subtelomere sequence that was shared between chromosomes in Amel_HAv3 (see Results) to help locating telomeres in this assembly using BLAST [40, 67]. Lastly, we queried Amel_HAv3 RepeatMasker output for shared interspersed repeats among the subtelomeric regions.

Satsuma v2 [44] was used to align Amel_4.5 against Amel_HAv3 using default settings. The alignments were used to characterize the sequences found to be shared between the assemblies, or unique to either of them. The Amel_4.5 sequence is not oriented with respect to the genetic map and we did therefore not perform in-depth assessments of synteny or reorientations between assemblies.

The mitochondrial genome sequence was recovered in a single contig. It was circularized and subject to two rounds of polishing. We then used BLAST to detect the order and orientation of the canonical set of coding sequences, rRNAs and tRNAs along the chromosome (NCBI accession NC_001566) [34, 42]. Because BLAST did not label all tRNAs in its default settings, we also used MITOS [41] together with MIFTI [75] to annotate the mitochondrial genome ab initio. The sequence was visualized in DNAPlotter [76]. In order to detect any major structural differences, we used MAFFT to align the whole mitochondrial sequence against the corresponding sequence in Amel_4.5.