# The black honey bee genome: insights on specific structural elements and a first step towards pan-genomes

Sonia Eynard, Christophe Kropp, Kamila Canale-Tabet, William Marande, Céline Vandecasteele, Céline Roques, Cécile Donnadieu, Quentin Boone, Bertrand Servin, Alain Vignal

► To cite this version:

Sonia Eynard, Christophe Kropp, Kamila Canale-Tabet, William Marande, Céline Vandecasteele, et al.. The black honey bee genome: insights on specific structural elements and a first step towards pan-genomes. 2024. hal-04473386

HAL Id: hal-04473386

https://hal.inrae.fr/hal-04473386

Preprint submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - NoDerivatives 4.0 International License
# The black honey bee genome: insights on specific structural elements and a first step towards pan-genomes

Sonia Eynard, Christophe Kropp, Kamila Canale-Tabet, William Marande, Céline Vandecasteele, Céline Roques, Cécile Donnadieu, Quentin Boone, Bertrand Servin, Alain Vignal

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sonia Eynard, Christophe Kropp, Kamila Canale-Tabet, William Marande, Céline Vandecasteele, et al.. The black honey bee genome: insights on specific structural elements and a first step towards pan-genomes. Genetics Selection Evolution, 2023, 56 (1), pp.51. 10.1101/2023.12.06.570386. hal-04473386v2

HAL Id: hal-04473386

https://hal.inrae.fr/hal-04473386v2

Submitted on 12 Jul 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
The black honey bee genome: insights on specific structural elements and a first step towards pan-genomes

3

4 Sonia E. Eynard¹, Christophe Klopp², Kamila Canale-Tabet¹, William Marande³, Céline

5 Vandecasteele⁴, Céline Roques⁴, Cécile Donnadieu⁴, Quentin Boone¹,², Bertrand Servin¹ and Alain

6 Vignal¹*

7 ¹GenPhySE, Université de Toulouse, INRAE, INPT, INP-ENVT, Castanet Tolosan, France

8 ²Sigenae, MIAT, INRAE, Castanet Tolosan, France

9 ³CNRGV, INRAE, Castanet Tolosan, France

10 ⁴INRAE, US 1426, GeT-PlaGe, Genotoul, Castanet-Tolosan, France

11

12 *Corresponding author

13

14 E-mail addresses:

15 SE: sonia.eynard@inrae.fr

16 CK: christophe.klopp@inrae.fr

17 KC-T: kamila.tabet@inrae.fr

18 WM: william.marande@inra.fr

19 CV: celine.vandecasteele@inra.fr

20 CR: celine.lopez-roques@inrae.fr

21 QB: quentin.boone@inrae.fr

22 CD: cecile.donnadieu@inrae.fr

23 BS: bertrand.servin@inrae.fr

24 AV: alain.vignal@inrae.fr

25
RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# The black honey bee genome: insights on specific structural elements and a first step towards pangenomes

Sonia E. Eynard<fnref n="1" />, Christophe Klopp<fnref n="2" />, Kamila Canale-Tabet<fnref n="1" />, William Marande<fnref n="3" />, Céline Vandecasteele<fnref n="4" />, Céline Roques<fnref n="4" />, Cécile Donnadieu<fnref n="4" />, Quentin Boone<fnref n="1" />,<fnref n="2" />, Bertrand Servin<fnref n="1" /> and Alain Vignal<fnref n="1" />*

**Abstract**

**Background** The honey bee reference genome, HAv3.1, was produced from a commercial line sample that was thought to have a largely dominant *Apis mellifera ligustica* genetic background. *Apis mellifera mellifera*, often referred to as the black bee, has a separate evolutionary history and is the original type in western and northern Europe. Growing interest in this subspecies for conservation and non-professional apicultural practices, together with the necessity of deciphering genome backgrounds in hybrids, triggered the necessity for a specific genome assembly. Moreover, having several high-quality genomes is becoming key for taking structural variations into account in pangenome analyses.

**Results** Pacific Bioscience technology long reads were produced from a single haploid black bee drone. Scaffolding contigs into chromosomes was done using a high-density genetic map. This allowed for re-estimation of the recombination rate, which was over-estimated in some previous studies due to mis-assemblies, which resulted in spurious inversions in the older reference genomes. The sequence continuity obtained was very high and the only limit towards continuous chromosome-wide sequences seemed to be due to tandem repeat arrays that were usually longer than 10 kb and that belonged to two main families, the 371 and 91 bp repeats, causing problems in the assembly process due to high internal sequence similarity. Our assembly was used together with the reference genome to genotype two structural variants by a pangenome graph approach with Graphyper2. Genotypes obtained were either correct or missing, when compared to an approach based on sequencing depth analysis, and genotyping rates were 89 and 76% for the two variants.

**Conclusions** Our new assembly for the *Apis mellifera mellifera* honey bee subspecies demonstrates the utility of multiple high-quality genomes for the genotyping of structural variants, with a test case on two insertions and deletions. It will therefore be an invaluable resource for future studies, for instance by including structural variants in GWAS. Having used a single haploid drone for sequencing allowed a refined analysis of very large tandem repeat arrays, raising the question of their function in the genome. High quality genome assemblies for multiple subspecies such as presented here, are crucial for emerging projects using pangenomes.

Alain Vignal
alain.vignal@inrae.fr

<footnote>*Correspondence:</footnote>

<footnote n="1">¹ GenPhySE, Université de Toulouse, INRAE, INPT, INP-ENVT, Castanet Tolosan, France</footnote>

<footnote n="2">² Sigenae, MIAT, INRAE, Castanet Tolosan, France</footnote>

<footnote n="3">³ CNRGV, INRAE, Castanet Tolosan, France</footnote>

<footnote n="4">⁴ INRAE, US 1426, GeT-PlaGe, Genotoul, Castanet-Tolosan, France</footnote>
26 **Abstract** (maximum 350 words).

27 **Background**

28 The actual honey bee reference genome, HAv3.1, was produced from a commercial line sample,
29 thought to have a largely dominant *Apis mellifera ligustica* genetic background. *Apis mellifera
30 mellifera*, often referred to as the black bee, has a separate evolutionary history and is the original
31 type in western and northern Europe. Growing interest in this subspecies for conservation and non-
32 professional apicultural practices, together with the necessity of deciphering genome backgrounds
33 in hybrids, triggered the necessity for a specific genome assembly. Moreover, having several high-
34 quality genomes is becoming key for taking structural variations into account in pan-genome
35 analyses.

36 **Results**

37 Pacific Bioscience technology long reads were produced from a single haploid black bee drone.
38 Scaffolding contigs into chromosomes was done using a high-density genetic map. This allowed for
39 a re-estimation of the honey recombination rate, over-estimated in some previous studies, due to
40 mis-assemblies resulting in spurious inversions in the older reference genomes. The sequence
41 continuity obtained is very high and the only limit towards continuous chromosome-wide sequences
42 seem to be due to tandem repeat arrays usually longer than 10 kb and belonging to two main
43 families, the 371 and 91 bp repeats, causing problems in the assembly process due to high internal
44 sequence similarity. Our assembly was used together with the reference genome, for genotyping
45 two structural variants by a pan-genome graph approach with Graphtyper2. Genotypes obtained
46 were either correct or missing, when compared to an approach based on sequencing depth analysis,
47 and genotyping rates were 89 and 76 % for the two variants respectively.

48 **Conclusions**
## Background

The honey bee *Apis mellifera* was originally found in Europe, Africa, and the Middle East, with the most eastern limit of its natural distribution situated in western Afghanistan, until a new subspecies was discovered in Kazakhstan [1]. The evolutionary origin of *Apis mellifera* is still unclear, with a possible origin in Eastern Africa or the Middle East, followed by colonization of Europe through different routes, leading to high genetic differentiation between geographically close populations or subspecies, namely *A. m. mellifera* (otherwise referred to as M-type) in western Europe on one side and *A. m. ligustica* from Italy or *A. m. carnica* (known as C-type) from eastern Europe on the other [2–5]. However, although *A. m. mellifera* is the original subspecies found in western Europe, it has become commonplace amongst breeders, in order to increase production or to facilitate the handling of colonies, to import other subspecies, mainly *A. m. ligustica* from Italy, *A. m. carnica* from Slovenia, and *A. m. caucasica* from Georgia. These are either to be bred as pure lines or as hybrids generated by artificial or directed insemination [6, 7]. As a consequence, these imported subspecies and hybrid lines mate naturally to local *A. m. mellifera* populations, threatening them and prompting the establishment of conservation programmes [8]. However, although it has been replaced in the majority of large professional beekeeper's facilities by imported honey bees, *A. m. mellifera* is still used by specialised breeders.

The honey bee reference genome, whose first version was obtained in 2006 [9], was updated twice: first in 2014 [10] and a second time in 2019, using long-read sequencing together with Hi-C chromatin interaction and BioNano Optical maps for a chromosome-scale assembly [11]. The sample used for this reference genome was from a commercial line (DH4) that was not precisely genetically defined, but is thought to be mainly of *A. m. ligustica* descent [9]. As a consequence, the genome of the genetically distinct *A. m. mellifera* may not be accurately represented and future pangenome approaches, which have been shown in other species to expand the number of genomic regions available for analysis [12, 13], would benefit from a high-quality assembly for this important subspecies.

To ensure a faithful representation of the *A. m. mellifera* subspecies genetic background, an individual from the black bee conservatory “Association Conservatoire de l’Abeille Noire Bretonne” in the island of Ouessant, France, was selected for sequencing. Ouessant is a very small island (15.5 km²), located 20 km off the coast of Brittany, in which a conservation population was set up starting in 1987. Any further imports of other honey bees are banned since 1991. Mitochondrial DNA analyses

have shown a low haplotype diversity and the presence of only the M-type in this population [14]. As expected from such a small population, microsatellite analysis has shown low diversity [15].

Until the latest update [11], the current honey bee genome sequence, Amel4.5 [10], suffered from imperfections, having numerous gaps in the assembly and possible sequence inversions. In order to construct a new *A. m. mellifera* genome assembly with improved continuity, we used the Pacific Biosciences long-read technology and produced all sequence reads from a single haploid drone to avoid assembly problems due to polymorphism. To order and orient our contigs along the chromosomes, we used published sequencing reads from drones originating from three colonies that had previously been used to map meiotic crossovers and non-crossovers in the honey bee [16], allowing also for the production of an updated genetic map and a re-estimation of the honey bee recombination rate.

Our analyses of the assembly allowed the detection of a major family of tandem repeats, running in some instances over more than 10 kb, and found at the ends of most sequence contigs. Our assembly allows for the first-time to perform detailed analyses of structural rearrangements, including at the population level, between the genomes of *A. m. ligustica* and other C-type honey bees that are used by the majority of beekeepers, and that of the M-type subspecies *A. m. mellifera* black bee.

## Methods

### Sampling, DNA extraction, and PacBio long-read sequencing

Candidate drones for sequencing were sampled at the larval or pupae stage from the black bee conservatory on the island of Ouessant, Brittany, France, and extractions were performed from several samples in order to select the best DNA quality in terms of molecular weight and quantity. Each sample was ground using a potter (see Additional file 1 Figure S1) and DNA extraction was performed using the QIAGEN Genomic-tips 100/G kit (Cat No./ID: 10243), following the tissue protocol extraction (see supplementary methods). DNA for sequencing was obtained from a single drone OUE7B (see Additional file 1 Figure S2). Library preparation and sequencing were performed at the GeT-PlaGe core facility, INRAE Toulouse, following the manufacturer's instructions for “Shared protocol-20 kb Template Preparation Using BluePippin Size Selection system (15 kb size Cutoff)”. At each step, DNA was quantified using the Qubit dsDNA HS Assay Kit (Life Technologies). DNA purity was tested using a nanodrop (Thermofisher) and size distribution and degradation was assessed using the Fragment analyzer (AATI) High Sensitivity Large Fragment 50 kb
49 Our new assembly for the *Apis mellifera mellifera* honey bee subspecies demonstrates the utility of
50 multiple high-quality genomes for the genotyping of structural variants, with a test case on two
51 insertions and deletions. It will therefore be an invaluable resource for future studies, for instance
52 including structural variants in GWAS. Having used a single haploid drone for sequencing allowed
53 a refined analysis of very large tandem repeat arrays, raising the question of their function in the
54 genome. High quality genome assemblies for multiple subspecies such as presented here, are crucial
55 for emerging projects using pan-genomes.

56

57

## Background

58 The honey bee *Apis mellifera* was originally found in Europe, Africa and the Middle East, with the
59 most eastern limit of its natural distribution situated in western Afghanistan until a new subspecies
60 was discovered in Kazakhstan [1]. The evolutionary origin of *Apis mellifera* is still unclear, with a
61 possible origin in Eastern Africa or the Middle East, followed by the colonization of Europe
62 through different routes, leading to high genetic differentiation between geographically close
63 populations or subspecies, namely *A. m. mellifera* (otherwise referred to as M-type) in western
64 Europe on one side and *A. m. ligustica* from Italy or *A. m. carnica* (known as C-type) from eastern
65 Europe on the other [2–5]. However, although *A. m. mellifera* is the original subspecies found in
66 western Europe, it has become commonplace amongst breeders, in order to increase production or
67 to facilitate the handling of colonies, to import other subspecies, mainly *A. m. ligustica* from Italy,
68 *A. m. carnica* from Slovenia and *A. m. caucasica* from Georgia, either to be bred as pure lines or as
69 hybrids generated by artificial or directed insemination [6,7]. As a consequence, these imported
70 subspecies and hybrid lines will mate naturally to local *A. m. mellifera* populations, threatening
71 them and prompting the establishment of conservation programmes [8]. However, although it has
Analysis Kit. Purification steps were performed using 0.45X AMPure PB beads (PacBio). Thirty µg of DNA was purified to perform three libraries. Using SMRTBell template Prep Kit 1.0 (PacBio), a DNA and end damage repair step was performed on 15 µg of unshared sample. Then blunt hairpin adapters were ligated to the libraries. The libraries were treated with an exonuclease cocktail to digest unligated DNA fragments. A size selection step using a 7 kb (library 1) or 9 kb (libraries 2 and 3) cutoff was performed on the BluePippin Size Selection system (Sage Science) with 0.75% agarose cassettes, Marker S1 high Pass 15-20 kb. Conditioned Sequencing Primer V2 was annealed to the size-selected SMRTbells. The annealed libraries were then bound to the P6-C4 polymerase using a ratio of polymerase to SMRTbell at 10:1. Then, after a magnetic bead-loading step (OCPW), SMRTbell libraries were sequenced on 36 SMRTcells on a RSII instrument from 0.05 to 0.2 nM, with a 360 min movie.

## Assembly into contigs and alignment to Amel4.5 for chromosome assignments

Raw reads were assembled with Canu 1.3 [17] using standard parameters and a first polishing of the assembly was done with quiver (version SMRT_Link v4.0.0) using standard parameters. The contigs obtained after the assembly step were aligned to the Amel4.5 reference genome using LAST v956 [18].

## Alignment of Illumina sequencing reads and SNP calling for crossing over analysis

All the Illumina paired-end sequences from Liu et al. [16] were downloaded from the NCBI SRA project SRP043350 (see Additional file 2 Table S1). The reads were aligned to the assembled contigs with BWA MEM v0.7.15 [19], duplicate reads were removed with Picard (v2.1.1; http://picard.sourceforge.net), and local realignment and base quality score recalibration (BQSR) was performed using GATKv3.7 [19]. SNPs were called in each drone independently with GATK HaplotypeCaller and consolidated into a single set of master sites, from which all individuals were genotyped with GATK GenotypeGVCFs (see scripts in supplementary material). Hard filtering of variants was performed with the following filters: FS > 60, MQ < 50 and SOR > 3; any SNP with missing genotypes were filtered out. Further quality controls were applied and, for each colony, SNPs falling into any of the following categories were discarded: (i) non-polymorphic SNPs in the colony, (ii) homozygous SNPs in the queen, (iii) heterozygous SNPs in drones, (iv) SNPs that appeared inconsistent with the observations in the two other colonies and (v) SNPs showing inconsistent allelic versions between queen and drone genotypes.

## Phasing and detection of recombination events

For each colony and informative SNP, genotyping results were used to define genotype vectors across all drones for the colony. Identical genotype vectors that followed one another within the same contig defined a segment with no observed crossing over in the drones of the colony and were grouped into bins. Not having access to grandparental genotypes, genotype phase between two successive bins within a contig was determined by finding which out of the two possible inverse vectors minimised the number of recombination events. Non-crossing-over gene conversion events, which can be misinterpreted as double recombination events and that occur usually on short DNA fragments, often shorter than a few kb [16], were removed to avoid inflating the size of the genetic map. Non-crossing-over gene conversion events were identified as: (i) bins of length shorter than 2 kb and occurring between two identical bins, or (ii) bins of length shorter than 2 kb for which the number of recombination events happening within this bin was higher than the number of recombination events needed to go from the bin before to the bin after it. Bins detected as non-crossing-over gene conversions were merged with their two identical surrounding bins. Both phasing and putative non-crossing over identification were performed iteratively from one bin to the next and independently for each colony. As a consequence, a set of phased vectors that minimised recombination events was obtained for each contig for each colony.

## Scaffolding contigs into chromosomes

Using the a priori assignment of contigs to chromosomes by alignment to Amel4.5 as a starting point, contigs were ordered and oriented iteratively in order to minimise the number of recombination events between the genotype vectors defined at their extremities. Contig scaffolding was first performed using the data for each colony separately and was thereafter confirmed using markers informative across all three colonies.

## Correction of the assembly with Illumina reads

Genomic DNA from the same individual as used for the PacBio sequencing was sequenced with an Illumina NovaSeq6000 instrument, producing over 28,000,000 reads (estimated raw sequencing depth = 37 X), NCBI SRA accession SRR15173860. These were aligned on the assembled genome with BWA MEM version 0.7.12-r1039 [20] using standard parameters. Variant detection was done with freebayes version 1.1.0 [21] and filtered to retain only those with a minimum quality score of 20 and '1/1' genotype or '0/1' with no read supporting the reference allele. Finally, corrections to the genome assembly
72 been replaced in the majority of large professional beekeeper's facilities by imported honey bees, *A. m. mellifera* is still used by dedicated breeders.

74 The honey bee reference genome, whose first version was obtained as soon as 2006 [9], was
75 updated twice: a first time in 2014 [10] and a second time in 2019, using long-read sequencing
76 together with Hi-C chromatin interaction and BioNano Optical maps for a chromosome-scale
77 assembly [11]. The sample used for this reference genome is from a commercial line (DH4), which
78 is not precisely genetically defined, but is thought to be mainly of *A. m. ligustica* descent [9]. As a
79 consequence, the genome of the genetically distinct *A. m. mellifera* may not be accurately
80 represented and future pangenome approaches, that were shown in other species to expand the
81 number of genomic regions available for analysis [12,13], would benefit from a high-quality
82 assembly for this important subspecies.

83 To ensure a faithful representation of the *A. m. mellifera* subspecies genetic background, an individual from the
84 black bee conservatory "Association Conservatoire de l'Abeille Noire Bretonne" in the island of Ouessant,
85 France was selected for sequencing. This very small island (15.5 km²) is located 20 km off the coast of Brittany,
86 the conservation population was set up starting in 1987, and further imports of other honey bees banned since
87 1991. Mitochondrial DNA analyses have shown a low haplotype diversity and the presence of only the M-type
88 in this population [14]. As expected from such a small population, microsatellite analysis has shown a low di-
89 versity [15].

90 Until the latest update [11], the current honey bee genome sequence, Amel4.5 [10] suffered from imperfections,
91 having numerous gaps in the assembly and possible sequence inversions. In order to construct a new *A. m. mel-
92 lifera* genome assembly with improved continuity, we used the Pacific Biosciences long-read technology and
93 produced all sequence reads from a single haploid drone to avoid assembly problems due to polymorphism. To
were done when alternative alleles were found in the VCF file using vcf-consensus from the vcftools package (version 0.1.12a) [22] with standard parameters.

## Comparison with Amel 4.5 and HAv3.1 assemblies

Estimation of recombination rate and positioning recombination events along the Amel4.5 and AMelMel1.1 assemblies was done using the same procedure as used for the de-novo assembly. GC content and sequence coverage for the queens' genotypes in AMelMel1.1 were measured in 0.5 Mb windows and the recombination rates were estimated using a script from Petit et al. [23] over 1 Mb windows. Completeness of the assemblies was estimated with BUSCO 3.0.2 [24] using OrthoDB v9.1 single-copy orthologs [25], from the Metazoa (n=978) and Hymenoptera (n=4415) BUSCO core set. Alignments of AMelMel1.1 to Amel4.5 and to HAv3.1 were done using LAST v956 [18]. Standard output psl files were produced to keep all alignments related to repeat elements, together with psl files from split alignments [18], corresponding to one-to-one alignments. Dotplot visualisation of alignments were produced with custom scripts, available at https://github.com/avignal5/Pacif icBee/tree/main. Inversions between the two genome assemblies were detected in the split alignment psl file. Liftovers of the HAv3.1 gtf and gff annotation to produce files with AMelMel1.1 annotation coordinates were done using CrossMap [26] and the chained alignment format output from the AMelMel1.1 to HAv3.1 LAST alignments.

## Analysis of repeat elements

Analysis of tandem repeats was done with Tandem Repeat Finder v4.09 (TRF) [27], setting the maximum period size to 2000 bp. The two major classes of repeat sizes, the 91 bp repeat and the 371 bp repeat, were analysed by aligning all repeats within a class size with MAFFT v7.313 [28]. Sequences reported by TRF from different parts of the genome start at different positions of the repeated element that was detected and to address this, the multifasta alignments produced by MAFFT were processed with a custom script to determine an identical arbitrary start point for all sequences before performing a second alignment with MAFFT. Phylogenetic trees were constructed using Jalview v2.11.2 [29] with the average distance option. Consensus sequences from all sequences selected within the groups defined based on the phylogenetic trees were used for a BLAST search in the AMelMel1.1 assembly and hits following one another at distances shorter than the repeat period size were grouped together. Finally, the previously described monomer consensus sequences, accession X57427.1 for AluI

and X89530.1 for *AvaI*, were used to detect their presence in the assembly by BLAST.

## Analysis of indels in populations

Indels were detected by aligning the two genomes HAv3.1 and AMelMAI1.1 to one another with minimap2 [30], followed by variant calling with SVIM-asm [31]. Two nuclear mitochondrial DNA (NUMT) were then selected for genotyping in a set of 80 haploid males, representing the three major European bee subspecies: *A. m. mellifera* (n=35), *A. m. ligustica* (n=30), and *A. m. caucasica* (n=15) (see Additional file 2 Table S2). All 80 samples were aligned to both assemblies, as described in Wragg et al. [6], and sequencing depth was estimated using SAMtools [32]. Individual genotypes in the sequencing data were determined for two selected indels by two methods. One method consisted of using GraphTyper2 [33] to detect breakpoints due to insertions, deletions, or inversions in the pangenome graph built with SVIM-asm using the two assemblies, HAv3.1 and AMelMel1.1. The other method consisted in using sequencing depths as an indication of presence or absence of indels. For a given indel and for each sample, the sequencing depth for the alignments on the genome in which the indel is present was calculated and compared to the sequencing depth of the sequences that flanked the indel on both sides. Normalisation was done by calculating the ratio between sequencing depth in the indel and in the flanking sequences. Determination of presence or absence of the indels was then done by K-means clustering with K=2. The vcf file produced by SVIM-asm was used to estimate the number and size distribution of indels that were larger than 40 bp. Bedtools v2.30.0 was used to check that indels did not overlap with contig boundaries in one or the other assembly and to search for overlaps between indels and exons in the GCF_003254395.2_ Amel_HAv3.1_genomic.gff annotation file that was downloaded from NCBI.

## Results

### PacBio long-read sequencing and assembly into contigs

All long-read sequence data came from a single haploid drone that was selected amongst several tested based on having the highest DNA concentration and a peak of DNA fragment length at 35 kb (see Additional file 1 Figure S2). A high proportion of reads exceeded 10 kb and a few reads were longer than 70 kb. Their size distribution is shown in Additional file 1 Figure S3 and S4. After assembly, a total of 200 contigs (gap-free sequence tracts) was obtained. The longest contig was 11.6 Mb and the N50 contig size was 5.1 Mb (see Additional file 2 Table S3 and see Additional file 1 Figure S5). These results are a major improvement in comparison to the 46 kb N50
94

95

96

97

98

99

100

101

102

order and orient our contigs along the chromosomes, we used published sequencing reads from drones originat-
ing from three colonies that had previously been used to map meiotic crossovers and non-crossovers in the
honey bee [16], allowing also for the production of an updated genetic map and a re-estimation of the honey
bee recombination rate.

Our analyses of the assembly allowed the detection of a major family of tandem repeats, running in some in-
stances over more than 10 kb and found at the ends of most sequence contigs. Our assembly allows for the
first-time to perform detailed analyses of structural rearrangements, including at the population level, between
the genomes of *A. m. ligustica* and other C-type honey bees used by the majority of beekeepers and that of the
M-type subspecies *A. m. mellifera* black bee.

103

Methods

104

*Sampling, DNA extraction and PacBio long-read sequencing.*

105

106

107

108

109

110

111

112

113

114

115

116

Candidate drones for sequencing were sampled at the larval or pupae stage from the black bee conservatory on the island of Ouessant, Brittany, France and extractions were performed from several samples, to select the best DNA quality in terms of molecular weight and quantity. Each sample was ground using a potter (see Additional file 1: Fig. S1) and DNA extraction performed using the QIAGEN Genomic-tips 100/G kit (Cat No./ID: 10243), following the tissue protocol extraction (see supplementary methods). DNA for sequencing was obtained from a single drone OUE7B (see Additional file 1: Fig. S2). Library preparation and sequencing were performed at the GeT-PlaGe core facility, INRAE Toulouse, following the manufacturer's instructions for “Shared protocol-20kb Template Preparation Using BluePippin Size Selection system (15kb size Cutoff)”. At each step, DNA was quantified using the Qubit dsDNA HS Assay Kit (Life Technologies). DNA purity was tested using the nanodrop (Thermofisher) and size distribution and degradation assessed using the Fragment analyzer (AATI) High Sensitivity Large Fragment 50kb Analysis Kit.
contig of Amel4.5 and quite similar to the N50 contig of
5.4 Mb observed in the HAv3.1 assembly [11]. Analysis
with BUSCO showed that overall, AMelMel1.1 had a
slightly larger gene content than both Amel 4.5 and the
most recently published reference assembly AmelHAv3.1
[11] (see Additional file 2 Table S4).

Chromosomal assignment and contig ordering
with crossing-over data

A priori chromosomal assignment of contigs was done
by alignment to the Amel4.5 assembly using LAST v956
[18]. Out of the 200 contigs, 110 aligned successfully.
Crossing-over data to confirm chromosome assign-
ment and the order of contigs along chromosomes was
obtained by using the reads from the sequencing of 43
drones from the three colonies that were initially used
to estimate recombination rate in honey bee [16]. Briefly,
this data set contains sequence data for three queens and
their drone offspring (15 to 13 depending on the colony).
Three of the drones of colony 1 were sequenced in dupli-
cate and were used for quality control of SNP calling.
Aligning these reads to our contigs allowed the detection
of 2,103,924 SNPs, on 176 contigs before quality control.
Out of these, approximately 64.5% were discarded due to
lack of polymorphism across the three colonies, 1% for
being homozygous in the queens, and 1% for being het-
erozygous in the drones. Furthermore, 0.2% of the SNPs
were discarded for being inconsistent between the three
drone replicates and 0.4% were discarded for having
allelic inconsistencies between queen and drones of the
same colony. After all the quality controls and for each
of the three colonies, 687,699, 698,123, and 672,728 reli-
able SNPs (approximately 32% of the initial SNPs), were
detected for the three colonies, on 114, 112, and 113
contigs respectively (see Additional file 1 Figure S6). In
total 120 contigs were at least partially informative across
the colonies, with 104 contigs informative in the three
colonies and 16 for only one or two. A total of 114,754
polymorphic SNPs was present in the 104 contigs that
were informative across all three colonies (see Additional
file 1 Figure S6). Genotype vectors for each SNP across
colony drones were then defined, allowing for detection
of within-contig cross-overs (see Additional file 1 Figure
S7). Genotype vectors from the ends of contigs were then
used to join contig ends together by finding for each con-
tig end, the best corresponding end from another contig
that had either the same genotype vector or a genotype
vector that presented a minimal number of cross-overs
(see Additional file 1 Figure S7). To minimize the number
of comparisons, the a priori chromosomal assignment by
alignment to Amel4.5 (see above) was used.

One hundred and two contigs out of the 110 with
chromosome assignment by sequence similarity to

Amel4.5 had SNP genotype data and were thus inform-
ative for cross-over detection. At least one cross-over
event, as evidenced by the presence of at least 2 geno-
type vector bins, could be detected within 86 of these
contigs, thus allowing for their orientation. The remain-
ing 16 contigs were oriented based on their alignment
to Amel4.5. All these contigs were small, except for
one contig on chromosome 7. For this contig, despite
its large size, close to 2.4 Mb, it was difficult to orien-
tate using the genetic map, as no cross-over could be
detected due to an unusually low number of SNPs and
a very low local recombination rate. Moreover, its ori-
entation could not be deduced from Amel4.5 or even
from the more recent assembly HAv3.1, as both pos-
sible orientations induced large inversions when com-
pared to these other two assemblies. Contigs assigned
to chromosomes by alignment only (8 contigs) or by
crossing-over data alone (16 contigs), were assigned to
their chromosomes, but at an unknown (unlocalised)
position. All remaining 72 contigs were considered
unplaced (see Additional file 1 Figure S6).

Tandem repeats at contig boundaries and orientation
of a large inversion on chromosome 7

With long-read data, sequence contigs are large but they
still don't cover the entire length of chromosomes, with
the exception of chromosome 16. When analysing the
contig ends, we found that almost all were composed of
tandem repeats arrays that were usually longer than the
read lengths, thus preventing assembly. To orientate the
large contig on chromosome 7, positioned as 5th in order
along the chromosome by the cross-over data, we took
advantage of the fact that the repeat elements detected by
TRF and that were present at both extremities of the con-
tig had different period sizes (258 and 1296 bp) and con-
sensus sequences. These were compared to the proximal
repeats of the 4th and the 6th contigs of chromosome
7. Interestingly, a tandem repeat element of 258 bp was
detected at the end of the 4th contig, and another one
of 1296 bp at the end of the 6th contig. Both had period
sizes that were identical to the extremities of the 5th con-
tig, suggesting the correct orientation of the 5th contig.
Correspondence between these contig ends was further
examined by pairwise alignment of the repeat sequences
with NCBI BLAST. Identity was 100% between the
sequences of identical period sizes, whereas no sig-
nificant similarity could be found between the others
(Fig. 1 and see Additional file 2 Table S5), thus confirm-
ing the orientation of the contig. Dotplots that compare
AMelMel1.1 and HAv3.1 are shown in Additional file 3
Figure S22 and suggest a very small number of discrepan-
cies, the major one residing on chromosome 7.
Purification steps were performed using 0.45X AMPure PB beads (PacBio). Thirty µg of DNA was purified to perform 3 libraries. Using SMRTBell template Prep Kit 1.0 (PacBio), a DNA and END damage repair step was performed on 15µg of unshared sample. Then blunt hairpin adapters were ligated to the libraries. The libraries were treated with an exonuclease cocktail to digest unligated DNA fragments. A size selection step using a 7kb (Library 1) or 9kb (libraries 2 and 3) cutoff was performed on the BluePippin Size Selection system (Sage Science) with the 0.75% agarose cassettes, Marker S1 high Pass 15-20kb. Conditioned Sequencing Primer V2 was annealed to the size-selected SMRTbells. The annealed libraries were then bound to the P6-C4 polymerase using a ratio of polymerase to SMRTbell at 10:1. Then after a magnetic bead-loading step (OCPW), SMRTbell libraries were sequenced on 36 SMRTcells on a RSII instrument from 0.05 to 0.2 nM with a 360 min movie.

*Assembly into contigs and alignment to Amel4.5 for chromosome assignments.*

Raw reads were assembled with Canu 1.3 [17] using standard parameters and a first polishing of the assembly was done with quiver (version SMRT Links v4.0.0) using standard parameters. The contigs obtained after the assembly step were aligned to the Amel4.5 reference genome using LAST v956 [18].

*Alignment of Illumina sequencing reads and SNP calling for crossing over analysis.*

All the Illumina paired-end sequences from Liu et al. (2015) [16] were downloaded from the NCBI SRA project SRP043350 (see Additional file 2: Table S1). The reads were aligned to the assembled contigs with BWA MEM v0.7.15 [19], duplicate reads removed with Picard (v2.1.1; http://picard.sourceforge.net), and local realignment and base quality score recalibration (BQSR) performed using GATKv3.7 [19]. SNPs were called in each drone independently with GATK HaplotypeCaller and consolidated into a single set of master sites, from which all individuals were
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Orientation of the AMelMel1.1 contig, presenting an inversion on chromosome 7 when compared to HAv3.1. The repeats present at the boundary between the contigs were used to orient the AMelMel1.1 contig on chromosome 7. Assemblies with one or the other orientation of the contig were self-aligned with LAST. Left: orientation from AMelMel1.0. Right: orientation from AMelMel1.1. For each pair of alignments, only the junction between contigs are shown: the two ends of the contig to orient, the end of the previous, and the start of next contigs. Results clearly show that the orientation in AMelMel1.1 is the correct one</figcaption></figure>

## Telomere and centromere consensus sequences

The presence of telomeres is an indication of completeness of the assembly. These were analysed by searching for the accepted TTAGG consensus sequence for Hymenoptera [34] in TRF analysis output, estimating their distance to the ends of chromosomes, and comparing the results to that of other 2–7 bp repeats, including non-TTAGG 5 bp repeats. Results (Fig. 2) showed that TTAGG were repeated with at least 842 copies when present at the extremities of chromosomes, whereas other interstitial TTAGG repeats had only 117 repeats or less (mean=21.3, median=16.7), a size distribution close to that of other pentanucleotide repeats (mean=24.2, median=14.4). See also Additional file 1 Figure S8 and Additional file 2 Tables S6, S7 and S8 for data on other STR motifs. In the AMelMel1.1 assembly, no TTAGG repeats were found on chromosomes 3, 7, 12, and 15, and they were found only at the beginning of chromosome 1. In contrast, in the HAv3.1 assembly, they could be found at both extremities of chromosome 1, but were absent

from chromosomes 5 and 11 [11]. An AATAT repeat was found at the beginning of chromosome 15 in our assembly.

The *AluI* and *AvaI* repetitive sequences, which were previously described as being respectively telomeric and centromeric [35], were localised on the AMelMel1.1 assembly by BLAST search and the number of copies per locus detected was counted (Fig. 2). The *AluI* repeat was found at the start of chromosomes 2 (6 repeats), 7 (3 repeats), 11 (46 repeats), and 12 (32 repeats). In addition, a single *AluI* element was found around position 8 Mb on chromosome 15, at more than 1.5 Mb from the distal end. Curiously, the *AluI* repeats found on chromosomes 2 and 11 were at the opposite end from the TTAGG sequences we detected (Fig. 2). The *AvaI* repeat was found as arrays at single loci on chromosomes 1, 2, 4, 9, 11, and 14. Only 4 copies in the array were found on chromosome 1, the other arrays having between 10 and more than 30 copies. The *AvaI* repeats were at the start of chromosomes 9 and 14, at the opposite end from the TTAGG repeats. On

(See figure on next page.)

**Fig. 2** Tandem repeats of period size 90–371 bp detected in the AMelMel1.1 assembly. The colour scale represents the period size of the repeat elements and the Y axis the total length of the repeat array. Vertical dotted lines represent the contig boundaries in the AMelMel1.1 assembly. The position of *AluI* and *AvaI* repeats are indicated with the number of repeats in parentheses. The figure shows clearly that most contigs are separated by tandem repeats of period size close to 371 bp, of length in the order of 10 kb or more. See also Additional file 1: Fig S11 for repeats of longer period size (1000–2000 bp). Although not represented on the graph (period size = 5), TTAGG telomere repeats are indicated with the number of repeats in parentheses, when present at a chromosome end
140 genotyped with GATK GenotypeGVCFs (see scripts in supplementary material). Any SNP with
141 missing genotypes were filtered out. Further quality controls were applied and for each colony,
142 SNPs falling into any of the following categories were discarded: i) non-polymorphic SNPs in the
143 colony, ii) homozygous SNPs in the queen, iii) heterozygous SNPs in drones, iv) SNPs that
144 appeared inconsistent with the observations in the two other colonies and v) SNPs showing
145 inconsistent allelic versions between queen and drone genotypes.

146 *Phasing and detection of recombination events.*

147 For each colony and informative SNP, genotyping results were used to define genotype vectors
148 across all drones for the colony. Identical genotype vectors following one another within a same
149 contig define a segment with no observed crossing over in the drones of the colony and were
150 grouped into bins. Not having access to grand-parental genotypes, genotype phase between two
151 successive bins within a contig was determined by finding which out of the two possible inverse
152 vectors minimised the number of recombination events. Non-crossing over gene conversion events,
153 which can be misinterpreted as double recombination events, occurring usually on short DNA
154 fragment often considered shorter than a few kb, [16] were removed to avoid inflating the size of
155 the genetic map. In our study, non-crossing over gene conversion events were identified as: i) bins
156 of length shorter than 2 kb, occurring between two identical bins, or ii) bins of length shorter than 2
157 kb for which the number of recombination events happening within this bin is higher than the
158 number of recombination events needed to go from the bin before to the bin after it. Bins detected
159 as non-crossing over gene conversions were merged with their two identical surrounding bins. Both
160 phasing and putative non-crossing over identification were performed iteratively from one bin to the
161 next and independently for each colony. As a consequence, a set of phased vectors minimising
162 recombination events was obtained for each contig in each colony.
Tandem repeats
Period size: 50-1000 bp. Nb repeats: 10-1000

<figure><img src="image_5.png" /><figcaption>**Fig. 2** (See legend on previous page.)</figcaption></figure>
163 Scaffolding contigs into chromosomes.

164 Using the a priori assignment of contigs to chromosomes by alignment to Amel4.5 as a starting
165 point, contigs were ordered and oriented iteratively in order to minimise the number of
166 recombination events between the genotype vectors defined at their extremities. The contig
167 scaffolding was first performed using the data for each colony separately and was thereafter
168 confirmed using markers informative across all three colonies.

169 *Correction of the assembly with Illumina reads.*

170 Genomic DNA from the same individual used for the PacBio sequencing was sequenced with an
171 Illumina NovaSeq6000 instrument, producing over 28 000 000 reads (estimated raw sequencing
172 depth = 37 X), NCBI SRA accession SRR15173860. These were aligned on the assembled genome
173 with BWA MEM version 0.7.12-r1039 [20] using standard parameters. Variant detection was done
174 with freebayes version 1.1.0 [21] and filtered to retain only those with a minimum quality score of
175 20 and '1/1' genotype or '0/1' with no read supporting the reference allele. Finally, corrections to the
176 genome assembly were done when alternative alleles were found in the VCF file using vcf-
177 consensus from the vcftools package (version 0.1.12a) [22] with standard parameters.

178 *Comparison with Amel 4.5 and HAv3.1 assemblies.*

179 Estimation of recombination rate and positioning recombination events along the Amel4.5 and
180 AMelMel1.1 assemblies was done following the same procedure as for the de-novo assembly. GC
181 content and sequence coverage for the queens' genotypes in AMelMel1.1 were measured in 0.5Mb
182 windows and the recombination rates were estimated using a script from Petit et al. (2017) [23] over
183 1Mb windows. Completeness of the assemblies was estimated with BUSCO 3.0.2 [24] using
184 OrthoDB v9.1 single-copy orthologs [25], from the Metazoa (n=978) and Hymenoptera (n=4415)
the other four chromosomes, they were at least at 1.8 Mb from a chromosome end (Fig. 2).

### Recombination pattern

Having used cross-over detection and a genetic map for contig scaffolding, we could estimate the total genetic map for AMelMel, which was approximately 50 Morgans long, giving an average recombination rate in the genome of 23 cM/Mb, close to the first estimates based on RAPD and microsatellite genetic maps [36–39] and to the most recent estimates based on SNPs [11, 40] (Table 1). However, although we used the same sequencing dataset as in Liu et al. [16], we found a drastic reduction in recombination rate between our genetic map and the one they initially published, which was 37 cM/Mb (Table 1). A big difference is that the latter was based on alignments of sequence reads to Amel4.5. When aligning our assembly to Amel4.5, we find an agreement on the chromosomal assignment of the contigs, but reveal many discrepancies in the orientation of large chromosome segments. At most breakpoint positions between the two assemblies, recombination hotspots were detected on Amel4.5 (Fig. 3 and see Additional file 4 Figure S23), suggesting these assembly errors were responsible for the overall higher recombination rate observed in Liu et al. [16]. This reduction from 37 to 23 cM/Mb is explained by these artefactual recombination hotspots detected in Amel4.5 at the breakpoint positions where the two assemblies disagree, that are absent in AMelMel1.1 (i.e. for chromosome 3 shown in Fig. 3 and see Additional file 4 Figure S23 for all the chromosomes).

### High conservation of tandem repeat sequences across chromosomes

We used TRF to further localise and analyse repeat arrays in the whole honey bee genome. Interestingly, two major period size classes for tandem repeats could be found: one in the size range of 91–93 bp, with a maximum number of 231 repeats, hereafter called the 91 bp repeat, and

the second in the size range of 367–371 bp, with a maximum number of 100 repeats, called the 371 bp repeat (see Additional file 1 Figure S9). The 91 repeats were found on all chromosomes, whereas the 371 bp repeats were on all chromosomes except chromosome 16 (see Additional file 1 Figure S10). Interestingly, very long repeats whose length is within the range of the sequence reads, were often found at the junction between two sequence contigs, confirming they could be responsible for the impossibility to sequence and assemble these regions properly (see Fig. 2, see Additional file 1 Figure S11).

We further investigated the nature of the 91 and 371 repeats by analysing their potential homogeneity in terms of sequence content. Summary statistics for the two classes showed very different distributions in terms of repeat copy numbers within tandem arrays (see Additional file 1 Figure S12 and see Additional file 2 Table S9). There was a total of 345 arrays of the 91 bp repeat in the genome and 131 arrays of the 371 bp repeats. However, these numbers dropped to 43 and 74 respectively when only considering tandem arrays of more than 10 repeats, suggesting that most of the 91 bp repeats had less than 10 elements (see Additional file 1 Figure S12). To investigate sequence homogeneity within each of the two repeat classes, we selected the repeat sequence defined by TRF for repeats that had strictly more than ten copies in tandem within an array. For the 91 bp repeat, we selected for $91 \le$ period size $\le$ 93 and for the 371 bp repeat $367 \le$ period size $\le$ 371, as suggested by the graph shown in Additional file 1 Figure S9. Then, for each repeat class, we performed a multi-sequence alignment with MAFFT, and produced an average distance tree with Jalview [29]. Results showed that out of the 74 sequences of the 371 bp repeat class, 72 were clearly grouped together, having high similarity (Fig. 4), whereas the 43 sequences of the 91 bp repeat class showed lower similarity. We therefore decided to subdivide the 91 bp repeat class into three groups of 20, 10, and 3

<table><caption>**Table 1** Literature comparison of Apis mellifera genetic maps</caption><thead><tr><th></th><th>Data</th><th>Physical size (Mb)</th><th>Genetic size (M)</th><th>CO/chromosome*</th><th>cM/Mb</th></tr></thead><tbody><tr><td>Hunt and Page [36]</td><td>Microsatellites</td><td>178</td><td>34.5</td><td>4.3</td><td>19.4</td></tr><tr><td>Solignac et al. [37]</td><td>Microsatellites</td><td>178</td><td>40.6</td><td>-</td><td>22.8</td></tr><tr><td>Solignac et al. [38]</td><td>Microsatellites</td><td>186</td><td>40</td><td>-</td><td>22.04</td></tr><tr><td>Beye et al. [39]</td><td>Microsatellites</td><td>238</td><td>45.5</td><td>5.7</td><td>19</td></tr><tr><td>Liu et al. [16]</td><td>SNP</td><td>220</td><td>81.4</td><td>5.1</td><td>37</td></tr><tr><td>Wallberg et al. [40]</td><td>SNP</td><td>229</td><td>59.5</td><td>-</td><td>26</td></tr><tr><td>Wallberg et al. [11]</td><td>SNP</td><td>219</td><td>47.3</td><td>-</td><td>21.6</td></tr><tr><td>AMelMel1.1</td><td>SNP</td><td>220</td><td>50</td><td>3.1</td><td>23</td></tr></tbody></table>

<footnote>*Number of cross-overs per chromosome</footnote>
BUSCO core set. Alignments of AMelMel1.1 to Amel4.5 and to HAv3.1 were done using LAST
v956 [18]. Standard output psl files were produced to keep all alignments related to repeat elements,
together with psl files from split alignments [18], corresponding to one-to-one alignments. Dotplot
visualisation of alignments were produced with custom scripts. Inversions between the two genome
assemblies were detected in the split alignment psl file. Liftovers of the HAv3.1 gtf and gff
annotation to produce files with AMelMel1.1 annotation coordinates were done using CrossMap
[26] and the chained alignment format output from the AMelMel1.1 to HAv3.1 LAST alignments.

Analysis of repeat elements.

Analysis of tandem repeats was done with Tandem Repeat Finder v4.09 (TRF) [27], setting the
maximum period size to 2000 bp. The two major classes of repeat sizes: the 91 bp repeat and the
371 bp repeat were analysed by aligning all repeats within a class size with MAFFT v7.313 [28].
Sequences reported by TRF from different parts of the genome start at different positions of the
repeated element detected and as a consequence, the multifasta alignments produced by MAFFT
were processed with a custom script, to determine an identical arbitrary start point for all sequences,
before performing a second alignment with MAFFT. Phylogenetic trees were calculated in Jalview
v2.11.2 [29] with the average distance option. Consensus sequences from all sequences selected
within the groups defined based on the phylogenetic trees were used for a BLAST search in the
AMelMel1.1 assembly and hits following one another at distances shorter than the repeat period
size were grouped together.

The previously described monomer consensus sequences: accession X57427.1 for *AluI* and
X89530.1 for *AvaI* were used to detect their presence in the assembly by BLAST.

Analysis of indels in populations.
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Comparison of Amel4.5 and AMelMel1.1 assemblies for chromosome 3. Abscissa: AMelMel. Ordinate: Amel4.5. AMelMel contig borders are represented with vertical dotted lines. In addition, for both Amel4.5 and AMelMel, the position and number of recombination events detected along the chromosome are represented for each interval flanked by informative markers in the meioses analyzed. Average SNP density and recombination rate are given for 1 Mb windows. Regions indicated in red on the Amel4.5 assembly represent recombination 'hotspots' regions, where the number of recombination events between two informative SNPs is higher than five. See supplementary data for the other chromosomes</figcaption></figure>

sequences, based on the average distance tree (Fig. 4). The remaining ten 91 bp repeat class sequences were singletons. A consensus sequence was made for each of the four group of sequences and used for a BLAST search in the AMelMel1.1 assembly. Homogeneity of the 371 bp consensus sequence was confirmed by the

detection of a very high number of hits of high similarity that covered the overall length of the queries (see Additional file 1 Figure S13). In contrast, for the three different consensus sequences that were used separately for the 91 bp repeat, alignment length and sequence similarities were lower, confirming that this repeat is
Indels were detected by aligning the two genomes HAv3.1 and AMelMA11.1 to one another with minimap2 [30], followed by variant calling with SVIM-asm [31]. Two nuclear mitochondrial DNA (NUMT) were then selected for genotyping in a set of 80 haploid males representing the three major European bee subspecies: *A. m. mellifera* (*n*=35), *A. m. ligustica* (*n*=30) and *A. m. caucasica* (*n*=15) (see Additional file 2: Table S2). All 80 samples were aligned to both assemblies as described in Wragg et al. (2022) [6] and sequencing depth was estimated using SAMtools [32]. Individual genotypes in the samples sequencing data was determined for the two indels by two methods. One method consisted in using GraphTyper2 [33], that will detect breakpoints due to insertions, deletions or inversions in the pangenome graph built with SVIM-asm using the two assemblies HAv3.1 and AMelMel11.1. The other method consisted in using sequencing depths as an indication of presence or absence of Indels. For a given Indel and for each sample, the sequencing depth for the alignments on the genome in which the Indel is present was calculated and compared to the sequencing depth of the sequences flanking the Indel on both sides. Normalisation was done by calculating the ratio between sequencing depth in the indel and in the flanking sequences.
Genotype presence or absence was then done by K-means clustering with K=2.

## Results

*PacBio long-read sequencing and assembly into contigs.*

All long-read sequence data comes from a single haploid drone selected amongst several tested, as having the highest DNA concentration and a peak of DNA fragment length at 35 kb (see Additional file 1: Fig. S2). A high proportion of reads exceeds 10 kb and a few reads are longer than 70 kb. Their size distribution is shown in Additional file 1: Fig. S3 and S4. After assembly, a total of 200 contigs (gap-free sequence tracts) was obtained. The longest contig is 11.6 Mb and the N50 contig size is 5.1 Mb (see Additional file 2: Table S3 and Additional file 1: Fig. S5). These results are a
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Phylogenetic trees for the tandem repeats of period size 91–93 and 367–371 bp. Only tandem repeats with ten or more elements, such as detected by Tandem Repeat Finder, were considered. Left: phylogenetic tree for the 74 sequences with a period size of 367–371 bp. Right: phylogenetic tree for the 43 sequences with a period size of 91–93 bp. The vertical red lines indicate the cut-off that was used to define the groups of sequence based on similarity</figcaption></figure>

defined more by its size, rather than by a specific repeat family based on sequence composition (see Additional file 1 Figure S13).

We then searched for the possible presence of the 371 and 91 bp repeats in other organisms. BLAST searches with each of the 371 bp repeat consensus sequences did not find any significant hit in the NCBI nucleic collection database. When searching with each of the 91 bp consensus repeats, four hits were found: three consensus sequences from repeat arrays from chromosome 11 and one consensus sequences from a repeat array from chromosome 12 showed sequence similarity to fragments of predicted lncRNAs LOC116185390, LOC105734921, LOC116415009 and LOC116185696, from unknown scaffolds of the genome assemblies of *Apis dorsata* and *Apis florea*. However, these lncRNAs are composed of two exons and span close to 1.5 kb in the genomes of *Apis dorsata* and *Apis florea*, suggesting that the 91 bp sequences correspond to only a portion (one out of two exons) of these lncRNAs. To investigate further, we performed BLAST searches with each of the consensus sequences directly on the refseq_genomes databases of *Apis cerana*, *Apis dorsata*, and *Apis florea* and a very high number of hits were found, suggesting that the 371 bp and 91 bp repeats were also present in these three

genomes, with an apparent slightly higher percent identity for the 91 bp repeat (see Additional file 1 Figure S14).

## Difference in the number of repeats of 5S ribosomal RNA genes

Genes that are repeated in tandem can often vary in numbers between individuals through unequal crossovers [41] and are, therefore, good candidates to study functional variation related to large rearrangements. A typical example of such genes is the 5S ribosomal RNA genes, whose copy number can vary greatly in the genome [42–44]. Alignment of a region from the AMelMel1.1 and HAv3.1 assemblies in a region on chromosome 3 that contained 5S ribosomal RNA genes, showed variation in the number of these genes between the two genomes (Fig. 5.). The period size of one of the repeat arrays of 5S genes was 357 bp, while that of the other was 373 bp. However, inclusion of this sequence in the multiple sequence analysis of the 371 bp repeat showed that these two sequences are different (see Additional file 1 Figure S15).

## Inversions and indels between AMelMel1.1 and HAv3.1

One-to-one split alignments produced by aligning AMelMel1.1 on HAv3.1 with LAST were used to detect
230 major improvement in comparison to the 46 kb N50 contig of Amel4.5 and quite similar to the N50
231 contig of 5.4 Mb observed in the HAv3.1 assembly [11]. Analysis with BUSCO showed that
232 overall, AMelMel1.1 had a slightly larger gene content than both Amel 4.5 and the most recently
233 published reference assembly AmelHAv3.1 [34] (see Additional file 2: Table S4).

234 Chromosomal assignation and ordering contigs with crossing-over data.

235 A priori chromosomal assignment of contigs was done by alignment to the Amel4.5 assembly using
236 LAST v956 [18]. Out of the 200 contigs, 110 aligned successfully. Crossing-over data for
237 confirming assignation and ordering contigs along chromosomes was obtained by using the reads
238 from the sequencing of 43 drones from three colonies, initially used to estimate recombination rate
239 in honey bee [16]. Briefly, this data set contains sequence for three queens and their drone offspring
240 (15 to 13 depending on the colony). Three of the drones of colony 1 are sequenced in duplicate and
241 are used for the quality control of SNP calling. Aligning these reads to our contigs allowed the
242 detection of 2,103,924 SNPs, on 176 contigs before any quality control. Out of these, approximately
243 64.5% were discarded due to an absence of polymorphism within the three colonies analysed, 1%
244 for being homozygous in the queens and 1% for being heterozygous in the drones. Furthermore,
245 0.2% of the SNPs were discarded for being inconsistent between the three drone replicates and
246 finally 0.4% were discarded for having allelic inconsistencies between queen and drones of the
247 same colony. After all the quality controls and for each of the three colonies, 687,699; 698,123 and
248 672,728 reliable SNPs (approximately 32% of the initial SNPs), were detected in each of the three
249 colonies on 114, 112 and 113 contigs respectively (see Additional file 1: Fig. S6). In total 120
250 contigs were at least partially informative across the colonies, with 104 contigs informative in the
251 three colonies and 16 for only one or two. A total of 114,754 polymorphic SNPs was present overall
252 in the 104 contigs informative across all three colonies (see Additional file 1: Fig. S6). Genotype
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Differences in copy numbers for 5S RNA ribosomal genes. Shown are two of the loci containing 5S RNA genes that are present at 15 kb distance on chromosome 3. Top: screenshot of the NCBI genome viewer for the region showing the annotation for the 5S RNA genes. Bottom: dotplot alignment of HAv3.1 (x-axis) and AMelMel1.1 (y-axis) in the region. The first group of genes in the bottom left contains seven genes in HAv3.1 and twenty in AMelMel1.1 on the forward strand. The second in the top right contains eleven genes in HAv3.1 and eight in AMelMel1.1 on the reverse strand. The red lines off diagonal show the sequence similarity between the two groups of genes and confirm that the two gene clusters are in reverse orientation</figcaption></figure>
vectors for each SNP across colony drones were then defined, allowing for within-contig crossing-over detection (see Additional file 1: Fig. S7). Genotype vectors from the ends of contigs were then used to join contig ends together by finding for each contig end, the best corresponding end from another contig having either the same genotype vector or a genotype vector presenting a minimal number of crossing-overs (see Additional file 1: Fig. S7). To minimize the number of comparisons, the *a priori* chromosomal assignment by alignment to Amel4.5 (see above) was used.

One hundred and two contigs out of the 110 with chromosome assignment by sequence similarity to Amel4.5, had SNP genotype data and were thus informative for crossing-over detection. At least one crossing-over event, as evidenced by the presence of at least 2 genotype vector bins, could be detected within 86 of these contigs, thus allowing for their orientation. The remaining 16 contigs were oriented based on the alignment to Amel4.5. All these contigs were small, except one contig on chromosome 7. Despite its large size, close to 2.4 Mb, it was indeed difficult to orientate using the genetic map, as no crossing-over could be detected due to an unusually low number of SNPs and a very low local recombination rate. Moreover, its orientation could not be deduced from Amel4.5 or even from the more recent assembly HAv3.1, as both possible orientations induced large inversions when compared to these other two assemblies. Contigs assigned to chromosomes by alignment only (8 contigs) or by crossing-over data alone (16 contigs), were assigned to their chromosomes, but at an unknown (unlocalised) position. All remaining 72 contigs were considered misplaced (see Additional file 1: Fig. S6).

*Tandem repeats at contig boundaries and Orientation of a large inversion on chromosome 7.*

With long read data, sequence contigs are large, but still don't cover the entire length of chromosomes, with the exception of chromosome 16. When analysing the contig ends, we found that almost all were composed of tandem repeats arrays usually longer than the read lengths, thus
inversions larger than 1000 bp between the two genomes. The largest inversion detected was on chromosome 7 and was larger than 1.6 Mb (see Additional file 3 Figure S22 and see Additional file 5 Figure S24). It should be noted, that a similar rearrangement on chromosome 7 was previously detected when comparing a genome assembly of an *A. m. ligustica* samples with the HAv3.1 reference [45]. Although close to one hundred other inversions could be detected, their visual inspection on dotplot graphs showed that 53 were within complex repeat patterns that were present at the junction between contigs, 32 within other complex repeat elements, and only 12 were in the middle of the high-quality sequence contigs in both assemblies, thus representing well supported inversions. Apart the large inversion on chromosome 7, the smallest was 1055 bp long and the largest 25,608 bp long (see Additional file 2 Table S10 and see Additional file 5 Figure S24). Interestingly, some inversions are in genes and can involve repeat elements found in both assemblies. In the example shown in Fig. 6, a local inverted duplicated region that was seen in the HAv3.1 assembly was absent in AMelMel1.1. This chromosomal segment contains a portion of the gene model *LOC113218640*, which has no direct annotation in the HAv3.1 assembly, but is described as coding for a *bric-a-brac 1-like* protein. *Bric-a-brac* was shown to be involved in body pigmentation in Drosophila [46]. Another interesting inversion is 11 kb long on chromosome 3, in an intron of *Rhomboid*, a gene involved in the formation of wing veins in Drosophila [47]. A more complex rearrangement involves a gene labelled as a probable nuclear hormone receptor, *HR38*, which is involved in synchronizing the reproductive activity in *Agrotis ipsilon* [48] and in the larvalpupal transition in *Leptinotarsa decemlineata* [49]. Other genes involved in the inversions described have various functions [50–55] and are reported in Additional file 2 Table S10.

The vcf file produced by SVIM-asm was used to study the size distribution of indels that were larger than 40 bp. In total, 9589 indels were detected, only 66 of which were larger than 1000 bp and six larger than 10 kb (see Additional file 1 Figure S16). The largest insertion of a segment that was present in AMelMel1.1 and absent in HAv3.1 that was detected was 31,640 bp long and the largest segment that was present in HAv3.1 and absent in AMelMel1.1 was 22,461 bp long. Two hundred and ninety-one indels overlapped totally or partially with exons, 289 of which were smaller than 725 bp (see Additional file 1 Figure S17). One insertion 17,759 bp long, specific to AMelMel, was located in exon 3 of the uncharacterised gene *LOC410644* and one 6,953 bp deletion on chromosome 10 removes exon 2 from the uncharacterised gene *LOC100577772* in AMelMel.

## Using both assemblies for analysis of two medium-size indels in honey bee subspecies

To demonstrate the utility of using two reference genomes for analysing structural variants, we studied two indels corresponding to nuclear mitochondrial DNA (NUMT), that were detected by using minimap2 [30] and SVIM-asm [31]. The first, NUMT_Chr2, is 745 bp long, has 92.7% identity over 99% of its length to HAv3.1 mitochondrial DNA, is present in the AMelMel1.1 assembly on chromosome 2 at positions 12,212,275–12,213,020, and is absent from the HAv3.1 assembly. The second, NUMT_Chr10, is 576 bp long, has 92.5% identity over 94% of its length to HAv3.1 mitochondrial DNA, is present in the HAv3.1 assembly on chromosome 10 at positions 670,675–671,251, and is absent in the AMelMel1.1 assembly. Presence and absence of these two NUMTs were tested in three honey bee subspecies: *A. m. mellifera* (n=35), *A. m. ligustica* (n=30) and *A. m. caucasia* (n=15), for which Illumina sequencing data was aligned to both reference genomes. Inspection of mean sequencing depth over all 80 samples in the regions of NUMT_Chr2 and NUMT_Chr10 indicated a decrease of the mean depth and an increase of its variance (see Additional file 1 Figure S18), suggesting the existence of a presence / absence polymorphism. When inspecting the sequencing depth per population, the *A. m. mellifera* samples showed a constant value over NUMT_Chr2 and a depth close to zero over NUMT_Chr10, whereas the *A. m. ligustica* show an inverse tendency (Fig. 7). The *A. m. caucasia* samples did not appear to have NUMT_Chr2 in their genomes, whereas a few may have NUMT_Chr10, as suggested by the incomplete drop of mean sequencing depth on HAv3.1 in the corresponding region (Fig. 7). To genotype our samples individually, we used two methods. The first was to estimate individual sequencing depth in the chromosomal region that delimits the NUMTs, by using AMelMel1.1 as reference genome for NUMT_Chr2 and HAv3.1 for NUMT_Chr10 (see methods). All 80 samples could thereafter be called unambiguously and assigned to one of two groups (presence or absence) by K-means clustering (see Additional file 1 Figure S19). The second method tested was to use GraphTyper2 [33], allowing the genotyping of structural variation using pangenome graphs. Our GraphTyper2 results showed that the calling of samples was incomplete, with a high proportion of no-calls, and that using individual bam files of alignments to one or to the other reference genome can greatly influence the call rate (see Additional file 2 Table S11). Indeed, for detection of variants with minimap2 and SVIM-asm, a reference genome must be specified and bam files of alignments to this specific reference genome must be used to perform individual genotyping. So, we first used HAv3.1 as reference, resulting
preventing assembly. To orientate the large contig on chromosome 7, positioned as 5th in order along the chromosome by the CO data, we took advantage of the fact that the repeat elements detected by TRF and present at both extremities of the contig have different period sizes (258 and 1296 bp) and consensus sequences. These were compared to the proximal repeats of the 4th and the 6th contigs of chromosome 7. Interestingly, a tandem repeat element of 258 bp was detected at the end of the 4th contig, and of 1296 bp at the end of the 6th contig, period sizes identical to the extremities of the 5th contig, suggesting the correct orientation of the 5th contig. Correspondence between these contig ends was further examined by pairwise alignment of the repeat sequences with NCBI BLAST. The Identity was 100 % between the sequences of identical period sizes, whereas no significant similarity could be found between the others (see Fig. 1 and Additional file 2: Table S9), thus confirming the orientation of the contig. Dotplots comparing AMelMel1.1 and HAv3.1 are shown in Additional file 3 and suggest a very small number of discrepancies, the major one residing on chromosome 7.

### *Telomeric and centromeric consensus sequences.*

The presence of telomeres is an indication of the completeness of the assembly. These were analysed by searching for the accepted TTAGG consensus sequence for Hymenoptera [35] in TRF analysis output, estimating their distance to the ends of chromosomes and comparing the results to that of other 2-7 bp repeats, including non-TTAGG 5 bp repeats. Results (Fig. 2) show that TTAGG are repeated with at least 842 copies when present at the extremities of chromosomes, whereas other interstitial TTAGG repeats have only 117 repeats or less (mean = 21.3, median = 16.7), a size distribution close to that of other pentanucleotide repeats (mean = 24.2, median = 14.4). See also Additional file 1: Fig. S8 and Additional file 2: Table S6, S7 and S8 for data on other STR motifs. In the AMelMel1.1 assembly, no TTAGG repeats were found on chromosomes 3, 7, 12 and 15 and
<figure><img src="image_9.png" /><figcaption>**Fig. 6** A 10 kb inverted duplication on chromosome 3 between HAv3.1 and AMelMel1.1. Bottom right: a dot plot representation of the alignment with LAST of AMelMel1.1 to HAv3.1 shows a 10 kb inversion on chromosome 3. Self-alignments of AMelMel1.1 (left) and HAv3.1 (top) show that the latter has an inverted repeated sequence in the region. The vertical yellow lines show the position of repeats that were previously detected and shown in the NCBI annotation (grey boxes) and that were also found in our LAST alignments. NCBI annotation of genes are in green</figcaption></figure>

in a genotyping call rate of 78.7% for NUMT_Chr2 and null for NUMT_Chr10, with no result in the output file from GraphTyper2 at all. To check if the reference genome could influence the results, we also performed the analysis by using AMelMel1.1 as reference and this

time the call rate was 85.0% for NUMT_Chr2, and 76.2% for NUMT_Chr10. When genotype calls were successfully obtained in both analyses, results were identical and also concordant with the analysis based on sequencing depth, showing that results were consistent when
299 were found only at the beginning of chromosome 1, whereas in the HAv3.1 assembly, they could be
300 found at both extremities of this chromosome, but were absent from chromosomes 5 and 11 [11].
301 An AATAT repeat was found at the beginning of chromosome 15 in our assembly.

302 The *AluI* and *AvaI* repetitive sequences, previously described as being respectively telomeric and
centromeric were localised on the AMelMel1.1 assembly by BLAST search and the number of
303 copies per locus detected was counted (Fig 2). The *AluI* repeat was found at the start of
304 chromosomes 2 (6 repeats), 7 (3 repeats), 11 (46 repeats) and 12 (32 repeats). In addition, a single
305 *AluI* element was found around position 8 Mb on chromosome 15, at more than 1.5 Mb from the
306 distal end. Curiously, the *AluI* repeats found on chromosomes 2 and 11 were at the opposite end
307 from the TTAGG sequences we detected (Fig. 2). The *AvaI* repeat was found as arrays at single loci
308 on chromosomes 1, 2, 4, 9, 11 and 14. Only 4 copies in the array were found on chromosome 1, the
309 other arrays having between 10 and more than 30 copies. The *AvaI* repeats are at the start of
310 chromosomes 9 and 14, at the opposite end from the TTAGG repeats. On the other four
311 chromosomes, they are at least at 1.8 Mb from a chromosome end (Fig. 2).
312

313 Recombination pattern

314 Having used crossing-over detection and a genetic map for contig scaffolding, we could estimate
315 the total genetic map for AMelMel, which is approximately 50 Morgans long, giving an average
316 recombination rate in the genome of 23 cM/Mb, close to the first estimations based on the
317 microsatellite genetic map and to the most recent ones based on SNPs (Table 1). However, although
318 we used the same sequencing dataset as in Liu et al. (2015) [16], we found a drastic reduction in
319 recombination rate between our genetic map and the one they initially published, which is 37
320 cM/Mb (Table 1). A great difference is that the latter is based on alignments of the sequence reads
321 on Amel4.5. When aligning our assembly with Amel4.5, we find an agreement on the chromosomal
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Insertions and deletions in *Apis mellifera* subspecies. Analysis of nuclear mitochondrial DNA (NUMT) insertions detected in only one assembly. Top: dotplot representation of LAST alignments between the two assemblies show a 745 bp variant present in AMelMel1.1 on chromosome 2 and absent in HAv3.1 (left) and a 576 bp variant present in HAv3.1 chromosome 10 and absent in AMelMel1.1 (right). For each variant, sequencing depths were evaluated for individual samples on the reference in which it is present. Bottom: mean sequencing depth per subspecies, with *A. m. caucasia* (15 samples) in green, *A. m. ligustica* (30 samples) in yellow and *A. m. mellifera* (35 samples) in black. Results suggest that most of the *A. m. mellifera* samples contain the insertion present in the AMelMel1.1 assembly on chromosome 2, as the sequencing depth remains constant throughout the region, and not the one present in the HAv3.1 assembly on chromosome 10, as indicated by a sequencing depth close to zero. Inversely, most of the *A. m. ligustica* samples contain the insertion present in the HAv3.1 assembly on chromosome 10 and not the one in the AMelMel1.1 assembly on chromosome 2. Most *A. m. caucasia* samples lack the insertion present in the AMelMel1.1 assembly and a few seem to have the insertion present in the HAv3.1 assembly</figcaption></figure>

genotyping was possible with Graphtyper2. Two samples were called as heterozygotes for NUMT_Chr2 when using AMelMel1.1 as reference and were counted as “no calls” because our samples were haploid. Low sequencing depth could have been a possible explanation for the absence of genotyping results with GraphTyper2 in some of the samples, but this does not seem to be the case, as all samples that failed genotyping had at least 8X average sequencing depth in the sequence flanking the NUMTS analysed, while successful genotyping could be obtained for samples having as little as 2X sequencing depth (see Additional file 1 Figure S20). Importantly, the individual genotyping results confirmed the overall impression that the presence or absence of the NUMT insertions are specific to the subspecies analysed, with most, if not all

samples having identical within-population genotypes, except for NUMT_Chr10 in *A. m. caucasia*, for which four out of eleven samples had a different allele. Interestingly, NUMT_Chr2 is present in all *A. m. mellifera* and in only two *A. m. ligustica* samples, and absent from all other samples, while NUMT_Chr10 is absent from *A. m. mellifera* samples and present in all but one *A. m. ligustica* samples and in four *A. m. caucasia* samples (Figs. 7, 8).

## Difference in sequencing depth according to the genome used for alignment

Overall alignment depths of sequence reads on both the HAv3.1 and AMelMel reference genome assemblies were determined for all 80 samples studied. The ratios between
assignment of the contigs, but reveal many discrepancies in the orientation of large chromosome segments. At most breakpoint positions between the two assemblies, recombination hotspots are detected on Amel4.5 (Fig. 3 and Additional file 4), suggesting these assembly errors were responsible for the overall higher recombination rate observed in Liu et al. (2015) [16]. This reduction from 37 cM/Mb to 23 cM/Mb is explained by these artefactual recombination hotspots detected on Amel4.5 at the breakpoint positions where the two assemblies disagree, that are absent in AMelMel1.1 (i.e. for chromosome 3 shown in Fig. 3 and Additional file 4 for all the chromosomes).

### *High conservation of tandem repeat sequences across chromosomes.*

We used TRF to further localise and analyse the repeat arrays in the whole honey bee genome.

Interestingly, two major period size classes for tandem repeats could be found: one in the size range of 91-93 bp, with a maximum number of 231 repeats, hereafter called the 91 bp repeat and the second in the size range of 367-371 bp, with a maximum number of 100 repeats, called the 371 bp repeat (see Additional file 1: Fig. S9). The 91 repeats are found on all chromosomes, whereas the 371 bp repeats are on all chromosomes except chromosome 16 (see Additional file 1: Fig. S10).

Interestingly, very long repeats whose length is within the range of the sequence reads, were often found at the junction between two sequence contigs, confirming they could be responsible for the impossibility to sequence and to assemble these regions properly (see Fig. 2, Additional file 1: Fig11).

We investigated further the nature of the 91 and 371 repeats by analysing their potential homogeneity in terms of sequence content. Summary statistics for the two classes show very different distributions in terms of repeat copy numbers within tandem arrays (see Additional file 1: Fig. S12 and Additional file 2: Table S9). There is a total of 345 arrays of the 91 bp repeat in the
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Comparing the indel variant calling between sequencing depth analysis and Graphyper2. Presence or absence of the NUMTs in the samples was evaluated by the pangenome graph approach with Graphyper2 (x-axis) and by estimating the sequencing depth at the position of the NUMTs on the genome in which they are present (y-axis). Sequencing depths were normalised by calculating the ratio between sequencing depth at the position of the NUMT sequence and that of the flanking sequence. Nine out of 80 samples (11%) could not be called for NUMT_Chr2 and 19 (24%) not for NUMT_Chr10. When alleles could be called by Graphyper2, results agreed with the data based on sequencing depth</figcaption></figure>

the alignment depths on each of the two genomes show that *A. m. caucasia* and *A. m. ligustica* samples had a slightly higher sequencing depth on HAv3.1 than on AMelMel1.1 and inversely, *A. m. mellifera* samples had a slightly higher sequencing depth on AMelMel1.1 than on HAv3.1 (see Additional file 1 Figure S21).

# Discussion

## AMelMel assembly quality and comparison to other honey bee assemblies

Although five chromosome level genome assemblies for *Apis mellifera* are available [56], ours has the originality of representing *Apis mellifera mellifera*. Indeed, this subspecies is genetically distinct from *Apis mellifera ligustica*, *Apis mellifera carnica*, and *Apis mellifera caucasia*

[6] which are represented by the four other assemblies. To ensure that the sample used is a thorough representative of the *A. m. mellifera* genetic background, our genome assembly relied on a single drone from a black honey bee conservation population that has been kept in isolation on a small island since 1991, thus preventing admixture. While this ensures at best that the AMelMel1.1 assembly is a true representative of the *A. m. mellifera* genetic ancestry, other populations from this subspecies in other geographic regions and also isolated for conservation purposes, may have diverged slightly from the sequenced individual. Future work should extend the repertoire of *A. m. mellifera* variation through long-read sequencing of samples from other origins. Another originality of our study is that the contigs we obtained were scaffolded
345 genome and 131 arrays of the 371 bp repeats. However, these numbers drop to 43 and 74
346 respectively when only considering tandem arrays of more than 10 repeats, suggesting that most of
347 the 91 bp repeats have less than 10 elements (see Additional file 1: Fig. S12). To investigate
348 sequence homogeneity within each of the two repeat classes, we selected the repeat sequence
349 defined by TRF for repeats having strictly more than ten copies in tandem within an array. For the
350 91 bp repeat, we selected for $91 \le$ period size $\le 93$ and for the 371 bp repeat $367 \le$ period size 371,
351 as suggested by the graph shown in Additional file 1: Fig. S9. Then, for each repeat class, we
352 performed a multi-sequence alignment with MAFFT, and produced an average distance tree with
353 Jalview. Results show that out of the 74 sequences of the 371 bp repeat class, 72 were clearly
354 grouped together, having high similarity (Fig. 4), whereas the 43 sequences of the 91 bp repeat class
355 showed lower similarity. We therefore decided to subdivide the 91 bp repeat class into three groups
356 of 20, 10 and 3 sequences, based on the average distance tree (Fig. 4). The remaining ten 91 bp
357 repeat class sequences were singletons. A consensus sequence was made for each of the four group
358 of sequences, and was used for a BLAST search in the AMelMel1.1 assembly. The homogeneity of
359 the 371 bp consensus sequence was confirmed by the detection of a very high number of hits of
360 high similarity covering the overall length of the queries (see Additional file 1: Fig. S13). On the
361 contrary, for the three different consensus sequences used separately for the 91 bp repeat, alignment
362 length and sequence similarities were lower, confirming that it to correspond more to a class size,
363 rather than a specific repeat family based also on sequence composition (see Additional file 1: Fig.
364 S13).
365

365 We then searched for the possible existence of the 371 and 91 bp repeats in other organisms.
366 BLAST searches with each of the 371 bp repeat consensus sequences did not allow to find any
367 significant hit in the NCBI nucleic collection database. When searching with each of the 91 bp
368 consensus repeats, four hits were found: three consensus sequences from repeat arrays from
into chromosomes using a genetic (recombination) map rather than the now more common HiC chromatin conformation and Bionano optical maps methods [57, 58]. Compared to the current HAv3.1 reference genome [11], our assembly is slightly longer (227 Mb versus 225 Mb) and was built from a lower number of contigs (200 versus 228) with very similar N50 contig values (5.1 Mb versus 5.4 Mb). However, the overall final coverage was slightly smaller (137 X Pac Bio and Illumina reads in AMelMel1.1 versus 192 X in HAv3.1). The BUSCO statistics were also very similar between AMelMel1.1 and HAv3.1 because the contig building was based on PacBio reads in both cases, with some correction using Illumina reads. In only one case did assembly of contigs into chromosomes using the recombination data fail to accurately order and orient a contig, a large contig on chromosome seven. Despite this limitation, we were able to orient this contig thanks to careful analysis of tandem repeat elements at its boundaries. Sequencing data for both HAv3.1 and our assembly, AmelMel1.1, are from a single haploid drone, which is a tremendous advantage for the resolution of regions that are largely composed of repeat elements. This was recently demonstrated in the human Telomere-to-Telomere project, for which a complete hydatidiform mole haploid cell line was used, which helped solve complex structures such as centromeres [59]. Our results show however that, although sequencing of repeat elements and especially of challenging tandem repeats seems to be resolved by the use of a single haploid sample and long reads, there are cases in which the total length of monotonous repeats is larger than the read lengths, preventing local assembly. As a result, for almost all contig boundaries investigated, long stretches of tandem repeats were found (Fig. 2). Interestingly, chromosome 16, which was obtained as a single contig, had no stretch of tandem repeats that exceeded 10 kb.

## Genetic maps and recombination rate in the honey bee

Having used genetic recombination data to scaffold our contigs, we were able to build a new recombination map and estimate 23 cM/Mbp for the overall recombination rate in the honey bee [16], which is of the same magnitude as the latest values from [11] and also congruent with prior values [38–40] (Table 1). It is interesting to note that the public sequencing dataset that we used, representing 43 drone genome offspring of three queens, gave a much higher estimate of 37 cM/Mb when it was previously used to generate genotyping data by alignment to the Amel4.5 reference genome [16]. On closer inspection, this higher overall recombination rate in Liu et al. [16] was due to very specific false recombination hotspots that appeared at contig junctions in Amel4.5, when at least one of them was inverted

compared to AMelMel1.1 (Fig. 3 and see Additional file 4 Figure S23). This illustrates the importance of the quality of the reference genome for such studies. Errors in local estimates of recombination rates when using a misassembled reference genome will in turn affect any analysis based on recombination maps or that include linkage disequilibrium.

## Tandem repeats and current limits for obtaining chromosome-wide contigs

We found a high occurrence of conserved tandem repeats in the honey bee genome, whose lengths and sequence conservation caused problems for scaffolding contigs into chromosomes, the ultimate goal being for each chromosome to be covered by a single contig. Indeed, long stretches of such repeats were found at the boundaries between contigs. Luckily, the only large contig in the assembly that could be placed on chromosome 7 but not oriented due to lack of sufficient genetic data, had different tandem repeats at each of its extremities, allowing us to decide on a correct orientation. However, other regions may still be problematic, the most striking example being the region between 1 and 3 Mb on chromosome 10. In this region, the contigs were small (<0.2 Mb) due to a high occurrence of tandem repeats, leading to difficulties in their ordering along the chromosome and their orientation. Moreover, these repeats appeared to mostly belong to the highly conserved 371 bp family, preventing their use for contig mapping. This portion of chromosome 10 has also been described as difficult to assemble in other studies [60].

## General chromosome structure: telomeres and centromeres

Cytogenetic studies based on fluorescent in situ hybridization of *AluI* and *AvaI* probes suggest that the honey bee genome is composed of one large metacentric and 15 acrocentric chromosomes [35]. This is to date still considered as the standard honey bee karyotype structure [11, 34]. However, other data could question this structure, for instance the suggested positions of the centromeres based on sequence characteristics of the HAv3.1 genome assembly, such as the (GC) content and the presence of *AluI* and *AvaI* repeats on chromosomes 7, 8, and 11 in Wallberg et al. [11].

Regarding telomeres, we were not able to identify the TTAGG consensus sequences on all 17 chromosome ends (two for the metacentric chromosome 1 and one for each of the other fifteen acrocentric chromosomes) where they were expected: none were detected on the right arm of chromosome 1 and on chromosomes 3, 12, and 15. Interestingly, some chromosomes also lacked TTAGG repeats in the HAv3.1 assembly, but these were
chromosome 11 and one consensus sequences from a repeat array from chromosome 12 showed
sequence similarity to fragments of predicted lncRNAs LOC116185390, LOC105734921,
LOC116415009, LOC116185696, from unknown scaffolds of the genome assemblies of Apis
dorsata and Apis florea. However, these lncRNAs are composed of two exons and span close to 1.5
kb in the genomes of Apis dorsata and Apis florea, suggesting that the 91 bp sequences correspond
to only a portion (one out of two exons) of these lncRNAs. To investigate further, we performed
BLAST searches with each of the consensus sequences directly on the refseq_genomes databases of
Apis cerana, Apis dorsata and Apis florea. A very high number of hits were found, suggesting the
371 bp and 91 bp repeats were also present in these three genomes, with an apparent slightly higher
percent identity for the 91 bp repeat (see Additional file 1: Fig. S14).

*Difference in the number of repeats of 5S ribosomal RNA genes.*

Genes that are repeated in tandem can often vary in numbers between individuals through unequal
crossing-over [36]. They are therefore good candidates to study functional variation related to large
rearrangements. A typical example of such genes is the 5S ribosomal RNA genes whose copy
number can vary greatly in the genome [37–39]. Alignment of a region from the AMelMel1.1 and
HAV3.1 assemblies in a region on chromosome 3 containing 5S ribosomal RNA genes, show a
variation in the number of these genes between the two genomes (Fig. 5.). The period size of one of
the repeat arrays of 5S genes is 357 bp, while that of the second is 373 bp. However, inclusion of
this sequence in the multiple sequence analysis of the 371 bp repeat shows that these two sequences
are different (see Additional file 1: Fig 15)

*Inversions between AMelMel1.1 and HAv3.1.*

One-to-one split alignments produced by aligning AMelMel1.1 on HAv3.1 with LAST were used to
detect inversions larger than 1000 bp between both genomes. The largest inversion detected is on
not the same as in ours (chromosomes 5 and 11). These discrepancies can be due to problems in the assembly of these repeat regions, either due to variations in sequence quality between the two datasets or to local variations in repeat content, rendering the assembly to be of varying difficulty due to biological reasons. It is interesting to note that in the older assemblies of the bee genome, based on the same DH4 strain as used for HAv3.1, extended analyses of telomeric and subtelomeric repeats showed that some chromosomes were easier to analyse than others and that no TTAGG repeats were identified for chromosome 5 [34]. Taken together, although the current sequencing data supports the consensus karyotype structure, we didn't find that the *AluI* repeat elements [35] could be considered as a marker of telomeres, because when such repeats were detected at the extremity of a chromosome, this was at the opposite end from the TTAGG repeats (see specifically chromosome 11 in Fig. 2).

The question of the exact position of the centromeres is a more complex one: the centromeres are expected to be at the middle of chromosome 1 and at the proximal end of each of the other chromosomes. The *Aval* repeat element, which is considered as a marker of centromeres [35], was not found on all chromosomes and even when found, the number of repeats in the array could be as small as four, such as on chromosome 1 (Fig. 2). With the exception of chromosome 11, for which an *Aval* repeat was found at position 5 Mb, the *Aval* elements, when they were detected on a chromosome, were found within 2.5 Mb of the chromosome ends, reflecting the results found for HAv3.1 [11]. However, although the positions of the *Aval* repeats was identical between the two assemblies, the number of repeat elements varied for each given position. Thus, for the moment, the exact position of the centromeres remains uncertain, but the presence of an *Aval* element remains a plausible indication, especially as these seem to be coincident with other specific characteristics, such as low (GC) content [40] or low levels of polymorphism and recombination rates [6]. If these characteristics are indicators of centromere positions, then chromosome 11 and perhaps also chromosome 7 should be considered sub-metacentric, although then TTAGG repeats would be expected at both of the extremities of these chromosomes, which is not the case in any of the studies to date. Further improvements in genome sequencing and assembly and in obtaining higher-resolution cytogenetic metaphase chromosome preparations will be necessary to elucidate this question.

## Comparing the genomes of two honey bee subspecies

The HAv3.1 assembly is based on a sample from the DH4 line, which is thought to be mainly of *A. m. ligustica*

descent [9]. The comparison with our *Apis mellifera melifera* AMelMel1.1 assembly allows for the detection of rearrangements that occurred between these two distinct genetic types that cannot be detected by short read sequencing. Short sequence fragments repeated in tandem, such as the 91 bp and 371 bp repeats described here, tend to vary in copy number through non-allelic homologous recombination (NAHR) or unequal cross-overs [61]. A rapid observation of the LAST alignment data between the two assemblies suggests that the 371 bp repeat element can vary greatly in copy number and the 91 bp element to a much lesser extent, although these preliminary observations will require more thorough analyses. To date, no obvious function was found for these elements, except for the fact that a BLAST search found that the 91 bp element shows similarity in sequence to one out of two exons of *Apis dorsata* and *Apis florea* lncRNAs, suggesting that these are incomplete and consequently not active in the repeat arrays. However, annotation of the lncRNAs in *Apis dorsata* and *Apis florea* is only based on alignment of short read RNA-seq data. More work is needed to confirm this finding concerning the 91 bp repeat and further comparisons with other bee genomes whose sequences are underway [62] will help understand these interesting genome elements.

The 5S ribosomal RNA genes are another interesting case of variation in gene number and studies in mouse and human have shown that this variation may be important for a balanced dosage of rRNA that can have implications in diseases [37, 38]. It would be interesting to determine whether the variations in 5S gene numbers observed here is a difference between the two honey bee subspecies investigated or if intra-population variation can be found.

After screening out rearrangements that could be due to errors associated with assembly problems, such as inversions of complete small contigs, thirteen inversions larger than 1 kb were detected between the two genomes. Of these, a large 1.6 Mb inversion on chromosome 7 is likely an error in HAv3.1, as it was also seen in the sequence of a closely related sample from the *Apis mellifera ligustica* subspecies [45]. Of the twelve remaining inversions, some involve genes that are present either at one of the breakpoints and that have inversions within their structure (usually introns) or whose structure remains intact but are in reverse orientation. Eventually, interesting functions that may explain some of the phenotypic differences found between the two subspecies represented by our dataset will be found (see Additional file 2 Table S10). Even when restricting to genes for which functions have been observed in insects, three genes stand out. One is *Bric-a-brac 1-like*, whose implication in body pigmentation in Drosophila [46], could
chromosome 7 and is larger than 1.6 Mb (see Additional file 3 and Additional file 5). It should be noted, that a similar rearrangement on chromosome 7 was previously detected when comparing a genome assembly of an *A. m. ligustica* samples with the HAv3.1 reference [40]. Although close to one hundred other inversions could be detected, their visual inspection on dotplot graphs show that 53 are within complex repeat patterns present at the junction between contigs, 32 within other complex repeat elements and only 12, are in the middle of the high-quality sequence contigs in both assemblies, thus representing well supported inversions. Apart the large inversion on chromosome 7, the smallest is 1055 bp long and the largest 25608 bp long (see Additional file 2: Table S10 and Additional file 5). Interestingly, some inversions will concern genes and can involve repeat elements differentially found in both assemblies. In the example shown in Fig. 6, a local inverted duplicated region seen in the HAv3.1 assembly, is absent in AMelMel1.1. This chromosomal segment contains a portion of the gene model LOC113218640, which has no direct annotation in the HAv3.1 assembly, but is described as coding for a *bric-a-brac 1-like* protein. *Bric-a-brac* was shown to be involved in body pigmentation in drosophila [41]. Another interesting inversion is 11 kb long on chromosome 3, in an intron of *Rhomboid*, a gene involved in the formation of wing veins in Drosophila [42]. A more complex rearrangement involves a gene labelled as a probable nuclear hormone receptor *HR38*, involved in synchronizing the reproductive activity in *Agrotis ipsilon* [43] and in the larval-pupal transition in *Leptinotarsa decemlineata* [44]. Other genes involved in the inversions described are reported in Additional file 2: Table S10.

*Using both assemblies for the analysis of two medium-size InDels in honey bee subspecies.*

To demonstrate the utility of using two reference genomes for analysing structural variants, we studied two indels corresponding to nuclear mitochondrial DNA (NUMT), that were detected by using minimap2 [30] and SVIM-asm [31]. The first one, NUMT_Chr2, is 745 bp long, has 92.7 %
be linked to our two reference genomes, representing light (yellow) and dark coloured honey bee subspecies. Another is *Rhomboid*, which has previously been shown to be involved in the formation of wing veins in Drosophila [47]. A third is the hormone receptor *HR38*, which has been shown to be involved in synchronisation of reproductive activity in the moth *Agrotis ipsilon* and in larval-pupal transition in the Colorado potato beetle *Lep tinotarsa decemlineata* [43, 44].

The number of insertions and deletions between the two assemblies detected by our approach was less than 10,000, only very few (66) were larger than 1000 bp, and fewer than 300 overlapped with annotated exons. The importance of these findings will have to be investigated in more detail, first by confirmation by long-read sequencing of additional samples, and also by improving the annotation of the assemblies, in which gene models may be imperfect, especially for non-coding genes.

## Perspectives for pangenomics and the test case of nuclear mitochondrial DNA segments

The limitations associated with the use of a single reference genome for an exhaustive analysis of all possible samples from a population has long been documented [63]. One of the most important limitations is the fact that sequence segments, often containing genes, can be specific to a limited number of individuals from a species or population. As a consequence, sequencing reads from such sequence segments that are present in a sample but absent in the reference, will at best be discarded or will map to wrong positions on the reference, thus creating false variants. Just by comparing two *A. mellifera* genome assemblies, we show here that large insertions and deletion can be found, 6 of which were larger than 10 kb. Although only a limited number overlapped with exons, their impact on phenotypes can be considerable through modifications of gene regulation, for instance by inserting, removing, or changing the position of enhancers and other regulatory elements. It is therefore essential that large structural variants be considered in association studies, and therefore be genotyped. Furthermore, loss of information through the use of a single reference genome can also concern SNPs, as in some regions, high densities of variants may lead to low mapping depths, due to a number of mismatches per read that exceed thresholds required for high mapping quality. We have indeed shown here that mapping reads from a given *A. mellifera* subspecies to a reference genome of the same subspecies yields higher sequencing depths than when a reference from a different subspecies is used. This suggests that some information is lost in the process. To address all these issues, the reference genome model can now be

extended to graph representations of genomes, in which alternate loci can be represented [64].

To test the potential of a pangenome approach for honey bees, we used the two reference genomes to genotype structural variants, in the form of two NUMTs that were present in either the HAv3.1 or in the AMelMel1.1 assembly. These variants were selected as a test case for presence or absence of a 745 bp fragment in the case of NUMT_Chr2 and a 576 bp fragment for NUMT_Chr10. Since the 80 population samples analysed here were all haploid, each of the NUMTs analysed should be either present or absent in a given individual, avoiding the uncertainties associated with heterozygote status and greatly simplifying the genotype calling. However, despite this, results showed that the Graphhtyper2 approach did not succeed in calling genotypes for all samples. Such results may be caused by the fact that Graphhtyper2 is biased towards one of these two genomes, although it used a vcf file that represented the alternative paths through the two genomes. This results potentially from the fact that (i) the population data needs to be used as bam files that resulted from a primary mapping to one or the other reference genome and (ii) the vcf file that represents the alternative paths through the genomes contains a standard representation of reference and alternate alleles. This reference-bias became really obvious when HAv3.1 was used as reference for the construction of the vcf and primary mapping of reads, NUMT_Chr10 could not be genotyped at all, while a higher genotype calling rate was obtained when AMelMel1.1 was used as the primary reference. Other methods for building graph genomes and aligning reads to them have been described, such as Minigraph-Cactus for graph construction [65] and Giraffe for the alignment of reads [66]. However, these remain biased, as a 'reference' assembly is chosen as an initial backbone, which is then augmented in turn with variation from the remaining assemblies [65]. Recently, new methods such as PanGenome Graph Builder (PGGB) [67] were introduced to address the reference-bias problem and could be tested. Although more exhaustive than single reference approaches, pangenomes are more complex to use in practice, due to specific problems they introduce, such as the definition of genome coordinates, which is essential when comparing different studies [66].

## Conclusions

In conclusion, we have presented a genome assembly for the honey bee *Apis mellifera mellifera* that is from a different subspecies than the current reference genome. One originality of the assembly process was the use of recombination data rather than optical maps
identity over 99 % of its length to HAv3.1 mitochondrial DNA, is present in the AMelMel1.1 assembly on chromosome 2 at positions 12,212,275 – 12,213,020 and absent from the HAv3.1 assembly. The second one, NUMT_Chr10, is 576 bp long, has 92.5 % identity over 94 % of its length to HAv3.1 mitochondrial DNA, is present in the HAv3.1 assembly on chromosome 10 at positions 670,675 – 671,251 and absent in the AMelMel1.1 assembly. The presence and absence of these two NUMTs were tested in three honey bee subspecies: *A. m. mellifera* (n=35), *A. m. ligustica* (n=30) and *A. m. caucasia* (n=15), for which Illumina sequencing data was aligned to both reference genomes. Inspection of mean sequencing depth over all 80 samples in the regions of NUMT_Chr2 and NUMT_Chr10 indicates a decrease of the mean depth and an increase of its variance (see Additional file 1: Fig. 16), suggesting the existence of a presence / absence polymorphism. When inspecting the sequencing depth per population, the *A. m. mellifera* samples show a constant value over NUMT_Chr2 and have a depth close to zero over NUMT_Chr10, whereas the *A. m. ligustica* show an inverse tendency (Fig. 7). The *A. m. caucasia* samples seem not to have NUMT_Chr2 in their genomes, whereas a few may have NUMT_Chr10, as although there is a drop of mean sequencing depth on HAv3.1 in the corresponding region, there is still some low coverage (Fig. 7). To genotype our samples individually, we used two methods. The first was to estimate individual sequencing depth in the chromosomal region delimiting the NUMTs, by using AMelMel1.1 as reference genome for NUMT_Chr2 and HAv3.1 for NUMT_Chr10 (see methods). All 80 samples could thereafter be called unambiguously assigned to one of two groups (presence or absence) by K-means clustering (see Additional file 1: Fig. 17). The second method tested was to use GraphTyper2 [33], allowing the genotyping of structural variation using pangenome graphs. Our GraphTyper2 results, showed that the calling of samples was incomplete, with a high proportion of no-calls, and that the fact of using individual bam files of alignments to one or to the other reference genome can greatly influence the call rate (see Additional file 2: Table S11). Indeed,
or HiC to scaffold contigs into chromosomes. We characterised for the first time long tandem repeats that are present in the genome and that were responsible for most sequence discontinuities and showed that these belong to two main repeats families that yet need to be further characterised and whose potential function in the genome remains to be investigated. Finally, we demonstrated the value of having two reference-quality genomes for the detection of structural variants, such as inversions and insertions-deletions and demonstrated the possibility of using a pangenome approach to genotype such variants in honey bee populations.

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12711-024-00917-3.

Additional file 1. Figure S1. Mortar and potter used for grinding individual drones. Figure S2. A fragment analyzer results for sample OUE7B. B fragment analyzer results for sample OUE8B. Figure S3. Distribution of read lengths per SMRTcell. Figure S4. Size distribution of sequencing reads from all 36 SMRTcells. Figure S5. Size distribution of the length of the contigs assembled with Canu. Figure S6. Alignment of contigs to Amel4.5 and detection of SNPs in the data from [16]. Figure S7. Ordering contigs by using crossing-over events. Figure S8. Telomeric repeats and other repeats of small period size. Figure S9. Period size and copy number of tandem repeats. Figure S10. Chromosome distribution of the 91 bp repeats and the 371 bp repeat. Figure S11. Tandem repeats of period size 1000-2000 bp detected in the AmelMel assembly. Figure S12. Distribution of the number of repeats in tandem arrays for the 91 bp and 371 bp repeats. Figure S13. Summaries from the BLAST hits on AmelMel1.1 obtained with the consensus sequences from the major group for the 371 bp repeats and the three major groups for the 91 bp repeats. Figure S14. BLAST searches with the 371 bp and the 91bp on *Apis* genomes. Figure S15. Phylogenetic trees for the tandem repeats of period size 367-374 bp. Figure S16. Sequencing depth within and around insertions and deletions in *Apis mellifera* subspecies. Figure S17. Sequencing depth and genotype deduction in NUMT_Chr2 and NUMT_Chr10. Figure S18. Sequencing depth and call rate with Graphtyper 2. Figure S19. Sequencing depth and genotype deduction in NUMT_Chr2 and NUMT_Chr10. Figure S20. Sequencing depth and call rate with Graphtyper 2. Figure S21. Overall alignment depth according to the reference used.

Additional file 2. Table S1. Illumina reads downloaded from NCBI SRA project SRP043350. Data from [16]. Table S2. Samples used for analysing the Indel variants. Table S3. Canu assembly statistics. Table S4. Statistics on the BUSCO genes expected to be found in metazoa and hymenoptera. Table S5. Repeats detected by Tandem Repeat Finder at the breakpoints of the 5th contig of chromosome 7. Table S6. Short tandem repeats within 1 kb of chromosome ends. Table S7. Short tandem repeats with > 400 copies. Table S8. Statistics on non-TTAGG Short Tandem Repeats with more than 10 repeat copies. Table S9. Summary statistics for the 91 and 371 bp repeats. Table S10. Inversions detected between AmelMel1.1 and HAv3.1. Table S11. Genotypes of NUMT_Chr2 and NUMT_Chr10 with Graphtyper and by examining local sequencing depth.

Additional file 3. Figure S22. Comparison of HAv3.1 and AmelMel1.1 genome assemblies. Each subpanel represents one of the 16 chromosomes of the honey bee.

Additional file 4. Figure S23. Comparison of Amel4.5 and AmelMel1.1 assemblies for each of the 16 chromosomes of the honey bee.

Additional file 5. Figure S24. Inversions larger than 1 kb detected between the AmelMel1.1 and HAv3.1 genome assemblies.

## Acknowledgements

We thank the Association Conservatoire de l'Abeille Noire Bretonne for the sample that was sequenced in this study.

## Author contributions

KC-T, WM and AV performed sampling and high molecular weight DNA extraction; CV, CR, CD performed the sequencing; CK performed assembly into contigs and contig quality checks; SEE, BS and AV did the chromosome-level assembly using recombination data; AV did the tandem repeats analyses; QB and AV did the NUMT detection and analysis. SEE and AV drafted the manuscript. All authors read and approved the final manuscript.

## Funding

This study was financially supported by the INRA Département de Génétique Animale (INRA Animal Genetics division) "PacificBee" grant. It was performed in collaboration with the GeT core facility, Toulouse, France (http://get.genotoul.fr), and was supported by France Génomique National infrastructure, funded as part of "Investissement d'avenir" program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09) and by the GET-PACBIO program (« Programme operationnel FEDER-FSE MIDI-PYRENEES ET GARONNE 2014-2020»).

## Availability of data and materials

The AmelMel1.1 assembly has been deposited on the NCBI under the accession number GCA_003314205. The reads of the 36 corresponding PACBIO_SMRT runs are in SRA under the accessions SRR9587836 to SRR9593684. Scripts and supplementary description of bioinformatic analyses are available in GitHub: https://github.com/avignal5/PacificBee/tree/main. Annotation files for the genome assembly are available at https://zenodo.org/records/10459067 for download.

## Declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

Received: 13 January 2024 Accepted: 4 June 2024

Published online: 28 June 2024

## References

1. Sheppard W, Meixner MD. *Apis mellifera pomonella*, a new honey bee subspecies from Central Asia. *Apidologie*. 2003;34:367–75.

2. Whitfield CW, Behura SK, Berlocher SH, Clark AG, Johnston JS, Sheppard WS, et al. Thrice out of Africa: ancient and recent expansions of the honey bee. *Apis mellifera*. *Science*. 2006;314:642–5.

3. Han F, Wallberg A, Webster MT. From where did the Western honeybee (*Apis mellifera*) originate? *Ecol Evol*. 2012;2:1949–57.

4. Cridland JM, Tsutsui ND, Ramírez SR. The complex demographic history and evolutionary origin of the western honey bee, *Apis mellifera*. *Genome Biol Evol*. 2017;9:457–72.

5. Dogantzis KA, Tiwari T, Conflitti IM, Dey A, Patch HM, Muli EM, et al. Thrice out of Asia and the adaptive radiation of the western honey bee. *Sci Adv*. 2021;7:2151.

6. Wragg D, Eynard SE, Basso B, Canale-Tabet K, Labarthe E, Bouchez O, et al. Complex population structure and haplotype patterns in the Western European honey bee from sequencing a large panel of haploid drones. *Mol Ecol Resour*. 2022;22:3068–86.

7. Pieplow JT, Brauße J, Van Praagh JP, Moritz RFA, Erler S. A scientific note on using large mixed sperm samples in instrumental insemination of honeybee queens. *Apidologie*. 2017;48:716–8.
439 for the detection of variants with minimap2 and SVIM-asm, a reference genome must be specified
440 and bam files of alignments to this specific reference genome must be used to perform the
441 individual genotyping. So, we first used HAv3.1 as reference and the results were a genotyping call
442 rate of 78.7 % for NUMT_Chr2 and null for NUMT_Chr10, as the line describing its potential
443 genotypes didn't appear in the output file from GraphTyper2 at all. To check if the reference
444 genome could influence the results, we also performed the analysis by using AMelMel1.1 as
445 reference and this time the call rate was 85.0 % for NUMT_Chr2, and 76.2 % for NUMT_Chr10.
446 When genotype calls were successfully obtained in both analyses, results were identical and were
447 also concordant with the analysis based on sequencing depth, showing that when genotyping was
448 possible with GraphTyper2, the results were consistent. Two samples were called as heterozygotes
449 for NUMT_Chr2, when using AMelMel1.1 as reference and were counted as "no calls", given our
450 samples were haploid. Low sequencing depth could have been a possible explanation for the
451 absence of genotyping results with GraphTyper2 in some of the samples, but this does not seem to
452 be the case, as all samples that failed genotyping had at least 8X average sequencing depth in the
453 sequence flanking the NUMTs analysed, whereas successful genotyping could be obtained for
454 samples having as little as 2X sequencing depth (see Additional file: Fig. 18). Substantially, the
455 individual genotyping results confirm the overall impression that the presence or absence of the
456 NUMT insertions are specific to the subspecies analysed, with most, if not all samples having
457 identical within-population genotypes, except for NUMT_Chr10 in A. m. caucasia, for which four
458 out of eleven samples have a different allele. Interestingly, NUMT_Chr2 is present in all A. m.
459 mellifera and only two A. m. ligustica samples, and absent from all other samples, whereas
460 NUMT_Chr10 is absent from A. m. mellifera samples and present in all but one A. m. ligustica
461 samples and four A. m. caucasia samples (Fig. 7, Fig. 8).
462

462
8. Pinto MA, Henriques D, Chávez-Galarza J, Kryger P, Garnery L, van der Zee R, et al. Genetic integrity of the Dark European honey bee (Apis mellifera mellifera) from protected populations: a genome-wide assessment using SNPs and mtDNA sequence data. J Apic Res. 2014;53:269–78.

9. Weinstock GM, Robinson GE, Worley KC, Hartfelder K, Zdobnov EM, Hartfelder K, et al. Insights into social insects from the genome of the honeybee *Apis mellifera*. Nature. 2006;443:931–49.

10. Elsik CG, Worley KC, Bennett AK, Beye M, Camara F, Childers CP, et al. Finding the missing honey bee genes: lessons learned from a genome upgrade. BMC Genomics. 2014;15:86.

11. Wallberg A, Bunikis I, Pettersson OV, Mosbech M-B, Childers AK, Evans JD, et al. A hybrid de novo genome assembly of the honeybee, *Apis mellifera*, with chromosome-length scaffolds. BMC Genomics. 2019;20:275.

12. Talenti A, Powell J, Hemmink JD, Cook EAJ, Wragg D, Jayaraman S, et al. A cattle graph genome incorporating global breed diversity. Nat Commun. 2022;13:910.

13. Wang T, Antonacci-Fulton L, Howe K, Lawson HA, Lucas JK, Phillippy AM, et al. The Human Pangenome Project: a global resource to map genomic diversity. Nature. 2022;604:437–46.

14. Garnery L, Franck P, Baudry E, Vautrin D, Cornuet J-M, Solignac M. Genetic diversity of the west European honey bee (Apis mellifera mellifera and *A. m. iberica*) I. Mitochondrial DNA. Genet Sel Evol. 1998;30:S31.

15. Garnery L, Franck P, Baudry E, Vautrin D, Cornuet J-M, Solignac M. Genetic diversity of the west European honey bee (Apis mellifera mellifera and *A. m. iberica*) II. Microsatellite loci. Genet Sel Evol. 1998;30:S49.

16. Liu H, Zhang X, Huang J, Chen J-Q, Tian D, Hurst LD, et al. Causes and consequences of crossing-over evidenced via a high-resolution recombination landscape of the honey bee. Genome Biol. 2015;16:15.

17. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 2017;27:722–36.

18. Frith MC, Kawaguchi R. Split-alignment of genomes finds orthologies more accurately. Genome Biol. 2015;16:106.

19. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.

20. Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010;26:589–95.

21. Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. arXiv. 2012. http://arxiv.org/abs/1207.3907.

22. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.

23. Petit M, Astruc J-M, Sarry J, Drouilhet L, Fabre S, Moreno CR, et al. Variation in recombination rate and its genetic determinism in sheep populations. Genetics. 2017;207:767–84.

24. Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO Applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35:543–8.

25. Zdobnov EM, Tegenfeldt F, Kuznetsov D, Waterhouse RM, Simão FA, Ioannidis P, et al. OrthoDB v9.1: cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs. Nucleic Acids Res. 2017;45:D744–9.

26. Zhao H, Sun Z, Wang J, Huang H, Kocher J-P, Wang L. CrossMap: a versatile tool for coordinate conversion between genome assemblies. Bioinformatics. 2014;30:1006–7.

27. Benson G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 1999;27:573–80.

28. Katoh K, Standley DM. MAFFT Multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.

29. Waterhouse AM, Procter JB, Martin DMA, Clamp M, Barton GJ. Jalview Version 2—a multiple sequence alignment editor and analysis workbench. Bioinformatics. 2009;25:1189–91.

30. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34:3094–100.

31. Heller D, Vingron M. SVIM-asm: structural variant detection from haploid and diploid genome assemblies. Bioinformatics. 2021;36:5519–21.

32. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

33. Eggertsson HP, Kristmundsdottir S, Beyter D, Jonsson H, Skuladottir A. Hardarson MT, et al. GraphTyper2 enables population-scale genotyping of structural variation using pangenome graphs. Nat Commun. 2019;10:5402.

34. Robertson HM, Gordon KHJ. Canonical TTAGG-repeat telomeres and telomerase in the honey bee, Apis mellifera. Genome Res. 2006;16:1345–51.

35. Beye M, Moritz RF. Characterization of honeybee (Apis mellifera L.) chromosomes using repetitive DNA probes and fluorescence in situ hybridization. J Hered. 1995;86:145–50.

36. Hunt GJ, Page RE. Linkage map of the honey bee, *Apis mellifera*, based on RAPD markers. Genetics. 1995;139:1371–82.

37. Solignac M, Vautrin D, Baudry E, Mougel F, Loiseau A, Cornuet J-M. A microsatellite-based linkage map of the honeybee, *Apis mellifera* L. Genetics. 2004;167:253–62.

38. Solignac M, Mougel F, Vautrin D, Monnerot M, Cornuet JM. A third-generation microsatellite-based linkage map of the honey bee, *Apis mellifera*, and its comparison with the sequence-based physical map. Genome Biol. 2007;8:R66.

39. Beye M, Gattermeier I, Hasselmann M, Gempe T, Schoiott M, Baines JF, et al. Exceptionally high levels of recombination across the honey bee genome. Genome Res. 2006;16:1339–44.

40. Wallberg A, Glémin S, Webster MT. Extreme recombination frequencies shape genome variation and evolution in the honeybee, *Apis mellifera*. PLoS Genet. 2015;11: e1005189.

41. Vignal A, London J, Rahuel C, Cartron JP. Promoter sequence and chromosomal organization of the genes encoding glycophorins A, B and E. Gene. 1990;95:289–93.

42. Gibbons JG, Branco AT, Godinho SA, Yu S, Lemos B. Concerted copy number variation balances ribosomal DNA dosage in human and mouse genomes. Proc Natl Acad Sci USA. 2015;112:2485–90.

43. Hall AN, Turner TN, Queitsch C. Thousands of high-quality sequencing samples fail to show meaningful correlation between 5S and 45S ribosomal DNA arrays in humans. Sci Rep. 2021;11:449.

44. Ding Q, Li R, Ren X, Chan L, Ho VWS, Xie D, et al. Genomic architecture of 5S rDNA cluster and its variations within and between species. BMC Genomics. 2022;23:238.

45. Cao L, Zhao X, Chen Y, Sun C. Chromosome-scale genome assembly of the high royal jelly-producing honeybees. Sci Data. 2021;8:302.

46. De Castro S, Peronnet F, Gilles J-F, Mouchel-Vielh E, Gibert J-M. *bric à brac* (*bab*), a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in *Drosophila melanogaster*. PLoS Genet. 2018;14: e1007573.

47. Sturtevant MA. The Drosophila rhomboid gene mediates the localized formation of wing veins and interacts genetically with components of the EGF-R signaling pathway. Genes Dev. 1993;7:961–73.

48. Gassias E, Durand N, Demondion E, Bourgeois T, Bozzolan F, Debernard S. The insect HR38 nuclear receptor, a member of the NR4A subfamily, is a synchronizer of reproductive activity in a moth. FEBS J. 2018;285:4019–40.

49. Shen C-H, Xu Q-Y, Mu L-L, Fu K-Y, Guo W-C, Li G-Q. Involvement of Leptinotarsa hormone receptor 38 in the larval-pupal transition. Gene. 2020;751: 144779.

50. Jiang H, Lv J. MicroRNA-301a-3p increases oxidative stress, inflammation and apoptosis in ox-LDL-induced HUVECs by targeting KLF7. Exp Ther Med. 2021;21:569.

51. Rotstein B, Molnar D, Adryan B, Llimargas M. Tramtrack is genetically upstream of genes controlling tracheal tube size in Drosophila. PLoS ONE. 2011;6: e28985.

52. Weiss F, Caruso V, De Rosa U, Beatino MF, Barbuti M, Nicoletti F, et al. The role of NMDA receptors in bipolar disorder: a systematic review. Bipolar Disord. 2023;25:624–36.

53. Durst MA, Ratia K, Lavie A. Identifying small molecule probes of ENTPD5 through high throughput screening. PLoS ONE. 2019;14: e0210305.

54. Boczonadi V, King MS, Smith AC, Olahova M, Bansagi B, Roos A, et al. Mitochondrial oxodicarboxylate carrier deficiency is associated with mitochondrial DNA depletion and spinal muscular atrophy-like disease. Genet Med. 2018;20:1224–35.

55. Janssens V, Goris J. Protein phosphatase 2A : a highly regulated family of serine/threonine phosphatases implicated in cell growth and signalling. Biochem J. 2001;53:417–39.

56. NCBI Genome. https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=7460. Accessed 05 April 2024.
Discussion 463

AMelMel assembly quality and comparison to other honey bee assemblies 464

Although five chromosome level genome assemblies for *Apis mellifera* are available [45] ours has the originality of representing *Apis mellifera mellifera*. Indeed, this subspecies is genetically distinct from *Apis mellifera ligustica*, *Apis mellifera carnica* and *Apis mellifera caucasia* [46] represented by the four other assemblies. Another originality of our study, is that the contigs we obtained were scaffolded into chromosomes using a genetic (recombination) map rather than the now more common HiC chromatin conformation and Bionano optical maps methods. Compared to the current HAv3.1 reference genome [11], our assembly is slightly longer (227 Mb versus 225 Mb), is built from a lower number of contigs (200 versus 228) with very similar N50 contig values (5.1 Mb versus 5.4 Mb). However, the overall final coverage was slightly smaller in our study (137X Pac Bio and Illumina reads versus 192X in HAv3.1). BUSCO statistics are also very similar due to the fact that contig building was based in both cases on PacBio reads with some correction using Illumina reads. Assembly of contigs into chromosomes using the recombination data failed to accurately order and orient in only one instance for a large contig on chromosome seven. Despite this limitation, the orientation of this contig was possible thanks to a careful analysis of tandem repeat elements at its boundaries. Sequencing data for both HAv3.1 and our assembly, AmelMel1.1, are from a single haploid drone, which is a tremendous advantage for the resolution of regions largely composed of repeat elements. This was recently demonstrated in the human Telomere-to-Telomere project, for which a complete hydatidiform mole haploid cell line was used, helping to solve complex structures such as centromeres [47]. Our results show however, that although the sequencing of repeat elements and especially of challenging tandem repeats seems resolved by the use of a single haploid sample and long reads, there are cases in which the total length of monotonous repeats is larger than the reads lengths, preventing local assembly. As a result, for
57. Putnam NH, O'Connell BL, Stites JC, Rice BJ, Blanchette M, Calef R, et al. Chromosome-scale shotgun assembly using an in vitro method for long-range linkage. *Genome Res.* 2016;26:342–50.

58. Shelton JM, Coleman MC, Herndon N, Lu N, Lam ET, Anantharaman T, et al. Tools and pipelines for BioNano data: molecule assembly pipeline and FASTA super scaffolding tool. *BMC Genomics.* 2015;16:734.

59. Nurk S, Koren S, Rhie A, Rautiainen M, Bzikadze AV, Mikheenko A, et al. The complete sequence of a human genome. *Science.* 2022;376:44–53.

60. Kaskinova M, Yunusbayev B, Altinbaev R, Raffiudin R, Carpenter MH, Kwon HW, et al. Improved *Apis mellifera* reference genome based on the alternative long-read-based assemblies. *G3 (Bethesda).* 2021;11:jkab223.

61. Hastings PJ, Lupski JR, Rosenberg SM, Ira G. Mechanisms of change in gene copy number. *Nat Rev Genet.* 2009;10:551–64.

62. Pennisi E. Sequencing all life captivates biologists. *Science.* 2017;355:894–5.

63. Rosenfeld JA, Mason CE, Smith TM. Limitations of the human reference genome for personalized genomics. *PLoS ONE.* 2012;7: e40294.

64. Church DM, Schneider VA, Steinberg KM, Schatz MC, Quinlan AR, Chin C-S, et al. Extending reference assembly models. *Genome Biol.* 2015;16:13.

65. Hickey G, Monlong J, Ebler J, Novak AM, Eizenga JM, Gao Y, et al. Pangenome graph construction from genome alignments with Minigraph-Cactus. *Nat Biotechnol.* 2024;42:663–73.

66. Sirén J, Monlong J, Chang X, Novak AM, Eizenga JM, Markello C, et al. Pangenomics enables genotyping of known structural variants in 5202 diverse genomes. *Science.* 2021;374:abg8871.

67. Garrison E, Guarracino A, Heumos S, Villani F, Bao Z, Tattini L, et al. Building pangenome graphs. *bioRxiv.* 2023. http://biorxiv.org/lookup/doi/https://doi.org/10.1101/2023.04.05.535718.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
487 almost all contig boundaries investigated, long stretches of tandem repeats were found (Fig. 2).

488 Interestingly, chromosome 16, which was obtained as a single contig, has no stretch of tandem repeats exceeding 10 kb.

490 Genetic maps and recombination rate in the honey bee

491 Having used genetic recombination data to scaffold our contigs, we could build a new recombination maps and give an estimation of 23 cM/Mbp for the overall recombination rate in the honey bee [16], which is of the same magnitude as the latest values from [34] and also congruent with prior values [48–50]. It is interesting to note, that the public sequencing dataset we used, representing 43 drone genome offspring of three queens, gave a much higher estimate of 37 cM/Mb when previously used for generating genotyping data by alignment on the Amel4.5 reference genome [16]. On closer inspection, this higher overall recombination rate in Liu et al. (2015) [16], is due to very specific false recombination hotspots that appear at contig junctions in Amel4.5, when at least one of them is inverted as compared to AMelMel1.1 (Fig. 3 and Additional file 4). This illustrates the importance of the quality of the reference genome for such studies. Errors in the local estimations of recombination rate when using mis-assembled reference genome will in turn affect any analysis based on recombination maps or including linkage disequilibrium.

503 Tandem repeats and the current limits for obtaining chromosome-wide contigs

504 We found a high occurrence of conserved tandem repeats in the honey bee genome, whose length and sequence conservation caused problems for scaffolding contigs into chromosomes, the ultimate goal being each chromosome covered by a single contig. Indeed, long stretches of such repeats were found at the boundaries between contigs. Luckily, the only large contig in the assembly, that could be placed on chromosome 7, but not oriented due to lack of sufficient genetic data, had different
tandem repeats at each of its extremities, allowing to decide on a correct orientation. However, other regions may still be problematic, the most striking example being the region between 1 and 3 Mb on chromosome 10. In this region, the contigs are small (< 0.2 Mb) due to a high occurrence of tandem repeats, leading to difficulties in their ordering along the chromosome and their orientation. Moreover, these repeats appear mostly to belong to the highly conserved 371 bp family, preventing their use for contig mapping. This portion of chromosome 10 has also been described as difficult to assemble in other studies [51].

General chromosome structure: telomeres, centromeres.

Cytogenetic studies based on fluorescent *in situ* hybridization of *AluI* and *AvaI* probes suggest that the honey bee genome is composed of one large metacentric and 15 acrocentric chromosomes [52]. This is to date still considered as the honey bee standard karyotype structure [34,35]. However, other data could question this structure, for instance the suggested positions of the centromeres based on sequence characteristics of the HAv3.1 genome assembly such as the (GC) % and the presence of *AluI* and *AvaI* repeats on chromosomes 7, 8 and 11 in Wallberg et al. (2019) [34].

Regarding telomeres, we were not able to identify the TTAGG consensus sequences on all the 17 chromosome ends (two for the metacentric chromosome 1 and one for each of the other fifteen acrocentric chromosomes) where they were expected: none were detected on the right arm of chromosome 1 and on chromosomes 3, 12 and 15. Interestingly, some chromosomes also lacked TTAGG repeats in the HAv3.1 assembly, but these were not the same (chromosomes 5 and 11).

These discrepancies can be due to problems in the assembly of these repeat regions, either due to variations in the sequence quality between the two datasets or to local variations in repeat content, rendering the assembly of varying difficulty due to biological reasons. It is interesting to note, that in the older assemblies of the bee genome, based on the same DH4 strain as HAv3.1, extended
analyses of telomeric and subtelomeric repeats showed that some chromosomes were easier to analyse than others and that no TTAGG repeats were identified for chromosome 5 [35]. Taken altogether, although the current sequencing data supports the actual consensus karyotype structure, we didn't find that the *AluI* repeat elements [52] could be considered as a marker of telomeres, as when such repeats were detected at the extremity of a chromosome, this was at the opposite end from the TTAGG repeats (see specifically chromosome 11 in Fig2).

The question of the exact position of the centromeres is a more complex one: the centromeres would be expected at the middle of chromosome 1 and at the proximal end of each of the other chromosomes. The *AvaI* repeat element, considered as a marker of the centromeres [52] was not found on all chromosomes and even when found, the number of repeats in the array could be as small as four, such as the repeat on chromosome 1 (Fig. 2). With the exception of chromosome 11 for which an *AvaI* repeat was found at the position 5 Mb, the *AvaI* elements, when detected on a chromosome, were found within 2.5 Mb of the chromosome ends, reflecting the results found on HAv3.1 [34]. However, although the positions of the *AvaI* repeats is identical between the two assemblies, the number of repeat elements vary for each given position. For the moment, the exact position of the centromeres remains uncertain, but the criteria of the eventual presence of an *AvaI* element remains a plausible indication, especially as these seem to be coincident with other specific characteristics, such as low (GC) content [50] or low levels of polymorphism and recombination rates [46]. If these characteristics are indicators of the centromere positions, then chromosome 11 and perhaps also chromosome 7 should be considered sub-metacentric, although in this case, TTAGG repeats would be expected at both of the extremities of these chromosomes, which is not the case in any of the studies to date. Further improvements in genome sequencing and assembly and in obtaining higher-resolution cytogenetic metaphase chromosome preparations will be necessary to elucidate this question.
556 *Comparing the genomes of two honey bee subspecies.*

557 The HAv3.1 assembly is based on a sample from the DH4 line, thought to be mainly of *A. m.*
558 *ligustica* descent [9]. The comparison with our *Apis mellifera mellifera* AMelMel1.1 assembly
559 allows for the detection of rearrangements occurring between these two distinct genetic types, that
560 can't be detected through short read sequencing.

561 Short sequence fragments repeated in tandem, such as the 91 bp and 371 bp repeats described here,
562 tend to vary in copy number through non-allelic homologous recombination (NAHR) or unequal
563 crossing-over [53]. A rapid observation of the LAST alignment data between the two assemblies
564 suggests that the 371 bp repeat element can vary greatly in copy number and the 91 bp element to a
565 much lesser extent, although these preliminary observations will require more thorough analyses.
566 No obvious function was found for these elements to date, except for the fact that a BLAST search
567 found that the 91 bp element shows similarity of sequence to one out of two exons of *Apis dorsata*
568 and *Apis florea* lncRNAs, suggesting these are incomplete and consequently not active in the repeat
569 arrays. However, the annotation of the lncRNAs in *Apis dorsata* and *Apis florea* is only based on
570 the alignment of short reads RNA-seq. More work is needed to confirm this finding concerning the
571 91 bp repeat and further comparisons with other bee genomes whose sequences are underway [54]
572 will help understand these interesting genome elements. The 5S ribosomal RNA genes are another
573 interesting case of variation in gene number and studies in mouse and human have shown that this
574 variation may be important for a balanced dosage of rRNA, that can have possible implications in
575 diseases [37,38]. It would be interesting to see if the variations of 5S gene numbers observed here is
576 a difference between the two honey bee subspecies investigated or if intra-population variation can
577 be found.
After screening out rearrangements that could be due to errors associated with assembly problems, such as inversions of complete small contigs, thirteen inversions larger than 1 kb were detected between the two genomes. Out of these, a large 1.6 Mb inversion on chromosome 7 is likely an error in HAv3.1, as it was also seen when sequencing a closely related sample from the *Apis mellifera ligustica* subspecies [40]. Out of the twelve remaining inversions, some involve genes, present either at one of the breakpoints, having inversions within their structure (usually introns) or whose structure remains intact, but are in reverse orientation. As usual, interesting functions that may explain some of the phenotypic differences found between the two subspecies represented by our dataset will be found (see Additional file 2: Table S10). Even when restricting to genes for which functions were observed in insects, three genes stand out. One is *Bric-a-brac 1-like*, whose implication in body pigmentation in *Drosophila* [41], could be linked to our two reference genomes representing light (yellow) and dark coloured honey bee subspecies. Another is *Rhomboid*, previously shown to be involved in the formation in wing veins in *Drosophila* [42]. A third is the hormone receptor *HR38*, shown to be involved in the synchronisation of reproductive activity in the moth *Agrotis ipsilon* and the larval-pupal transition in the Colorado potato beetle *Leptinotarsa decemlineata* [43,44].

*Nuclear mitochondrial DNA segments and perspectives for pangenomics.*

To test the utility of having two reference genomes for genotyping structural variants, we tried genotyping two NUMTs, present in one or the other HAv3.1 and AMelMel1.1 assembly, with Graphyper2. Results show that Graphyper2 could not call genotypes for all samples. In the first instance, this is surprising, given the fact that this test of presence or absence of a 745 bp fragment in the case of NUMT_Chr2 and a 576 bp one for NUMT_Chr10 is done on haploid samples, simplifying the problem, as each of the NUMTs should be either present or absent in each
individual tested. This may be caused by the fact that Graphyper2 extracts reads that were
previously mapped to the structural variant regions on a linear reference genome, thus possibly
introducing a bias. It is however surprising, that when HAv3.1 was used as reference for the
primary mapping of reads, NUMT_Chr10 could not be genotyped at all. This reference-bias could
be overcome by using more recent methods in which the reads for the genomes to genotype are
mapped directly on the pan-genome graph, although such methods are more complex to use in
practice, due to problems such as the definition of genome coordinates [55].

608

609

Conclusions

In conclusion, we present here a genome assembly for the honey bee *Apis mellifera* that is from a
different subspecies than the current reference genome. One originality of the assembly process was
to use recombination data rather than optical maps or HiC for scaffolding contigs into
chromosomes. We characterise for the first time long tandem repeats that are present in the genome
and are responsible for most sequence discontinuities and show that these belong to two main
repeats families yet to be further characterised and whose potential function in the genome remains
to be investigated. Finally, we show the interest of having two reference-quality genomes for the
detection of structural variants, such as inversions and insertions-deletions and demonstrate the
possibility of using a pan-genome approach for genotyping such variants in honey bee populations.

619

620

Declarations

621 Ethics approval and consent to participate

622 Not applicable
623 Consent for publication

624 Not applicable

625 Availability of data and materials

626 The AMelMel1.1 assembly has been deposited on the NCBI under the accession number

627 GCA_003314205. The reads of the 36 corresponding PACBIO_SMRT runs are in SRA under the

628 accessions SRR9587836 to SRR9593684. Scripts and supplementary description of bioinformatic

629 analyses are available in GitHub: https://github.com/avignal5/PacificBee/tree/main.

630 Competing interests

631 The authors declare that they have no competing interests

632 Funding

633 This study was financially supported by the INRA Département de Génétique Animale (INRA

634 Animal Genetics division) “PacificBee” grant. It was performed in collaboration with the GeT core

635 facility, Toulouse, France (http://get.genotoul.fr), and was supported by France Génomique

636 National infrastructure, funded as part of “Investissement d'avenir” program managed by Agence

637 Nationale pour la Recherche (contract ANR-10-INBS-09) and by the GET-PACBIO program («

638 Programme operationnel FEDER-FSE MIDI-PYRENEES ET GARONNE 2014-2020 »).

639 Authors' contributions

640 KC-T, WM and AV performed sampling and high molecular weight DNA extraction; CV, CR, CD

641 performed the sequencing; CK performed assembly into contigs and contig quality checks; SEE, BS

642 and AV did the chromosome-level assembly using recombination data; AV did the tandem repeats

643 analyses; QB and AV did the NUMT detection and analysis. SEE and AV drafted the manuscript.

644 All authors read and approved the final manuscript.
## References

647 1. Sheppard W, Meixner MD. Apis mellifera pomonella, a new honey bee subspecies from Central Asia.
648 Apidologie. 2003;34:367-75.

649 2. Whitfield CW, Behura SK, Berlocher SH, Clark AG, Johnston JS, Sheppard WS, et al. Thrice out of
650 Africa: ancient and recent expansions of the honey bee, Apis mellifera. Science. 2006;314:642-5.

651 3. Han F, Wallberg A, Webster MT. From where did the Western honeybee (Apis mellifera) originate? Ecol
652 Evol. 2012;2:1949-57.

653 4. Cridland JM, Tsutsui ND, Ramírez SR. The Complex Demographic History and Evolutionary Origin of
654 the Western Honey Bee, Apis Mellifera. Genome Biology and Evolution. 2017;9:457-72.

655 5. Dogantzis KA, Tiwari T, Conflitti IM, Dey A, Patch HM, Muli EM, et al. Thrice out of Asia and the
656 adaptive radiation of the western honey bee. Sci Adv. 2021;7:eabj2151.

657 6. Wragg D, Eynard SE, Basso B, Canale-Tabet K, Labarthe E, Bouchez O, et al. Complex population
658 structure and haplotype patterns in the Western European honey bee from sequencing a large panel of
659 haploid drones. Molecular Ecology Resources. 2022;22:3068-86.

660 7. Pieplow JT, Brauße J, Praagh JP, Moritz RFA, Erler S. A scientific note on using large mixed sperm
661 samples in instrumental insemination of honeybee queens. 2017;1-3.

662 8. Pinto MA, Henriques D, Chávez-Galarza J, Kryger P, Garnery L, van der Zee R, et al. Genetic integrity of
663 the Dark European honey bee (Apis mellifera mellifera) from protected populations: a genome-wide
664 assessment using SNPs and mtDNA sequence data. Journal of Apicultural Research. 2014;53:269-78.

665 9. Weinstock GM, Robinson GE, Worley KC, Hartfelder K, Zdobnov EM, Hartfelder K, et al. Insights into
666 social insects from the genome of the honeybee Apis mellifera. Nature. 2006;443:931-49.

667 10. Elsik CG, Worley KC, Bennett AK, Beye M, Camara F, Childers CP, et al. Finding the missing honey
668 bee genes: lessons learned from a genome upgrade. BMC Genomics. 2014;15:86.

669 11. Wallberg A, Bunikis I, Vinnere Pettersson O, Mosbech M-B, Childers AK, Evans JD, et al. A hybrid de
670 novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds. bioRxiv.
671 2018;1-37.

672 12. Talenti A, Powell J, Hemmink JD, Cook EAJ, Wragg D, Jayaraman S, et al. A cattle graph genome
673 incorporating global breed diversity. Nat Commun. 2022;13:910.

674 13. Wang T, Antonacci-Fulton L, Howe K, Lawson HA, Lucas JK, Phillippy AM, et al. The Human
675 Pangenome Project: a global resource to map genomic diversity. Nature. 2022;604:437-46.

676 14. Garnery L, Franck P, Baudry E, Vautrin D, Cornuet J-M, Solignac M. Genetic diversity of the west
677 European honey bee (Apis mellifera mellifera and A. m. iberica) I. Mitochondrial DNA. Genetics Selection
678 Evolution. 1998;30:S31.
679 15. Garnery L, Franck P, Baudry E, Vautrin D, Cornuet J-M, Solignac M. Genetic diversity of the west European honey bee (Apis mellifera mellifera and A. m. iberica) II. Microsatellite loci. Genetics Selection Evolution. 1998;30:S49.

682 16. Liu H, Zhang X, Huang J, Chen J-Q, Tian D, Hurst LD, et al. Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee. Genome Biol. 2015;16:15.

684 17. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive *k*-mer weighting and repeat separation. Genome Res. 2017;27:722-36.

686 18. Frith MC, Kawaguchi R. Split-alignment of genomes finds orthologies more accurately. Genome Biol. 2015;16:1-17.

688 19. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297-303.

691 20. Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010;26:589-95.

693 21. Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. arXiv:12073907 [q-bio] [Internet]. 2012 [cited 2019 Apr 3]; Available from: http://arxiv.org/abs/1207.3907

695 22. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156-8.

697 23. Petit M, Astruc J-M, Sarry J, Drouilhet L, Fabre S, Moreno C, et al. Variation in Recombination Rate and Its Genetic Determinism in Sheep Populations. Genetics. 2017;genetics.300123.2017.

699 24. Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO Applications from Quality Assessments to Gene Prediction and Phylogenomics. Molecular Biology and Evolution. 2018;35:543-8.

702 25. Zdobnov EM, Tegenfeldt F, Kuznetsov D, Waterhouse RM, Simão FA, Ioannidis P, et al. OrthoDB v9.1: cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs. Nucleic Acids Research. 2017;45:D744-9.

705 26. Zhao H, Sun Z, Wang J, Huang H, Kocher J-P, Wang L. CrossMap: a versatile tool for coordinate conversion between genome assemblies. Bioinformatics. 2014;30:1006-7.

707 27. Benson G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Research. 1999;27:573-80.

709 28. Katoh K, Standley DM. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. Molecular Biology and Evolution. 2013;30:772-80.

711 29. Waterhouse AM, Procter JB, Martin DMA, Clamp M, Barton GJ. Jalview Version 2--a multiple sequence alignment editor and analysis workbench. Bioinformatics. 2009;25:1189-91.

713 30. Li H. Minimap2: pairwise alignment for nucleotide sequences. Birol I, editor. Bioinformatics. 2018;34:3094-100.
715 31. Heller D, Vingron M. SVIM-asm: structural variant detection from haploid and diploid genome assemblies. Robinson P, editor. Bioinformatics. 2021;36:5519-21.

717 32. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078-9.

719 33. Eggertsson HP, Kristmundsdottir S, Beyter D, Jonsson H, Skuladottir A, Hardarson MT, et al. GraphTyper2 enables population-scale genotyping of structural variation using pangenome graphs. Nat Commun. 2019;10:5402.

722 34. Wallberg A, Bunikis I, Pettersson OV, Mosbech M-B, Childers AK, Evans JD, et al. A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds. BMC Genomics. 2019;20:275.

725 35. Robertson HM, Gordon KHJ. Canonical TTAGG-repeat telomeres and telomerase in the honey bee, Apis mellifera. Genome Res. 2006;16:1345-51.

727 36. Vignal A, London J, Rahuel C, Cartron JP. Promoter sequence and chromosomal organization of the genes encoding glycophorins A, B and E. Gene. 1990;95:289-93.

729 37. Gibbons JG, Branco AT, Godinho SA, Yu S, Lemos B. Concerted copy number variation balances ribosomal DNA dosage in human and mouse genomes. Proc Natl Acad Sci USA. 2015;112:2485-90.

731 38. Hall AN, Turner TN, Queitsch C. Thousands of high-quality sequencing samples fail to show meaningful correlation between 5S and 45S ribosomal DNA arrays in humans. Sci Rep. 2021;11:449.

733 39. Ding Q, Li R, Ren X, Chan L, Ho VWS, Xie D, et al. Genomic architecture of 5S rDNA cluster and its variations within and between species. BMC Genomics. 2022;23:238.

735 40. Cao L, Zhao X, Chen Y, Sun C. Chromosome-scale genome assembly of the high royal jelly-producing honeybees. Sci Data. 2021;8:302.

737 41. De Castro S, Peronnet F, Gilles J-F, Mouchel-Vielh E, Gibert J-M. bric à brac (bab), a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in Drosophila melanogaster. Kopp A, editor. PLoS Genet. 2018;14:e1007573.

740 42. Sturtevant MA. The Drosophila rhomboid gene mediates the localized formation of wing veins and interacts genetically with components of the EGF-R signaling pathway. Genes & Development. 1993;7:961-73.

743 43. Gassias E, Durand N, Demondion E, Bourgeois T, Bozzolan F, Debernard S. The insect HR38 nuclear receptor, a member of the NR4A subfamily, is a synchronizer of reproductive activity in a moth. FEBS J. 2018;285:4019-40.

746 44. Shen C-H, Xu Q-Y, Mu L-L, Fu K-Y, Guo W-C, Li G-Q. Involvement of Leptinotarsa hormone receptor 38 in the larval-pupal transition. Gene. 2020;751:144779.

748 45. NCBI Genome [Internet]. [cited 2023 Aug 28]. Available from: https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=7460

750 46. Wragg D, Eynard SE, Basso B, Canale-Tabet K, Labarthe E, Bouchez O, et al. Complex population
751 structure and haplotype patterns in Western Europe honey bee from sequencing a large panel of haploid

752 drones [Internet]. Genetics; 2021 Sep. Available from:

753 http://biorxiv.org/lookup/doi/10.1101/2021.09.20.460798

754 47. Nurk S, Koren S, Rhie A, Rautiainen M, Bzikadze AV, Mikheenko A, et al. The complete sequence of a

755 human genome. Science. 2022;376:44–53.

756 48. Beye M, Gattermeier I, Hasselmann M, Gempe T, Schioett M, Baines JF, et al. Exceptionally high levels

757 of recombination across the honey bee genome. Genome Res. 2006;16:1339–44.

758 49. Solignac M, Mougel F, Vautrin D, Monnerot M, Cornuet JM. A third-generation microsatellite-based

759 linkage map of the honey bee, Apis mellifera, and its comparison with the sequence-based physical map.

760 Genome Biol. 2007;8:R66.

761 50. Wallberg A, Glémin S, Webster MT. Extreme Recombination Frequencies Shape Genome Variation and

762 Evolution in the Honeybee, Apis mellifera. PLoS Genet. 2015;11:e1005189.

763 51. Kaskinova M, Yunusbayev B, Altinbaev R, Raffiudin R, Carpenter MH, Kwon HW, et al. Improved *Apis*

764 *mellifera* reference genome based on the alternative long-read-based assemblies. G3

765 Genes|Genomes|Genetics. 2021;11:jkab223.

766 52. Beye M, Moritz RF. Characterization of honeybee (Apis mellifera L.) chromosomes using repetitive

767 DNA probes and fluorescence in situ hybridization. J Hered. 1995;86:145–50.

768 53. Hastings PJ, Lupski JR, Rosenberg SM, Ira G. Mechanisms of change in gene copy number. Nat Rev

769 Genet. 2009;10:551–64.

770 54. Pennisi E. Sequencing all life captivates biologists. Science. 2017;355:894–5.

771 55. Sirén J, Monlong J, Chang X, Novak AM, Eizenga JM, Markello C, et al. Pangenomics enables

772 genotyping of known structural variants in 5202 diverse genomes. Science. 2021;374:abg8871.

773

774

775
776 **Figures**

777 **Figure 1. Orientation of the AMelMel1.1 contig presenting an inversion on chromosome 7**
778 when compared to HAv3.1.

779 The repeats present at the boundary between the contigs were used to orient the AMelMel1.1 contig
780 on chromosome 7. Assemblies with one or the other orientation of the contig were self-aligned with
781 LAST. Left: orientation from AMelMel1.0 and right, orientation from AMelMel1.1. For each pair
782 of alignments, only the junction between contigs are shown: the two ends of the contig to orient, the
783 end of the previous and the start of next contigs. Results clearly show the orientation in
784 AMelMel1.1 is the correct one.

785 **Figure 2. Comparison of Amel4.5 and AMelMel1.1 assemblies for chromosome 3.**

786 Abscissa: AMelMel, ordinate: Amel4.5. AMelMel contig borders are represented with vertical
787 dotted lines. Additionally, for both Amel4.5 and AMelMel, the position and number of
788 recombination events detected along the chromosome are represented for each interval flanked by
789 informative markers in the meiosis analyzed. Average SNP density and recombination rate are
790 given for 1Mb windows. Regions indicated in red on the Amel4.5 assembly represent
791 recombination ‘hotspots’ regions where number of recombination events between two informative
792 SNPs is higher than five. See supplementary data for the other chromosomes.

793 **Figure 3. Tandem repeats of period size 90-371 bp detected in the AMelMel1.1 assembly.** The
794 colour scale represents the period size of the repeat elements and the Y axis the total length of the
795 repeat array. Vertical dotted lines represent the contig boundaries in the AMelMel1.1 assembly. The
796 position of *AluI* and *AvaI* repeats are indicated with the number of repeats in parentheses. The
797 figure shows clearly, that most contigs are separated by tandem repeats of period size close to 371
798 bp, of length in the order of 10 kb or more. See also Supplementary file 1: Fig S11 for repeats of
799 longer period size. (1000-2000 bp). Although not represented on the graph (period size = 5),
800 TTAGG telomere are indicated with the number of repeats in parentheses, when present at a
801 chromosome end.

802 **Figure 4. Phylogenetic trees for the tandem repeats of period size 91-93 and 367-371 bp.**

803 Only tandem repeats with ten or more elements, such as detected by Tandem Repeat Finder, were
804 considered. Left: phylogenetic tree for the 74 sequences with a period size of 367-371 bp; right:
805 phylogenetic tree for the 43 sequences with a period size of 91-93 bp. The vertical red lines indicate
806 the cut-off that was used to define the groups of sequence based on similarity.

807 **Figure 5. Differences in copy numbers for 5S RNA ribosomal genes.** Two of the loci containing
808 5S RNA genes, present at 15 kb distance on chromosome 3 are shown. Top: screenshot of the NCBI
809 genome viewer for the region showing the annotation for the 5S RNA genes. Bottom: dotplot
810 alignment of HAv3.1 (x-axis) and AMelMel1.1 (y-axis) in the region. The first group of genes in
811 the bottom left contains seven genes in HAv3.1 and twenty in AMelMel1.1 on the forward strand.
812 The second in the top right contains eleven genes in HAv3.1 and eight in AMelMel1.1 on the
813 reverse strand. The red lines off diagonal show the sequence similarity between the two groups of
814 genes and indicate the two gene clusters are in reverse orientation.

815 **Figure 6. A 10 kb inverted duplication on chromosome 3 between HAv3.1 and AMelMel1.1.**

816 Bottom right: a dot plot representation of the alignment with LAST of AMelMel1.1 to HAv3.1
817 show a 10 kb inversion on chromosome 3. Self-alignments of AMelMel1.1 (left) and HAv3.1 (top)
818 show that the latter has an inverted repeated sequence in the region. The vertical yellow lines show
819 the position of repeats that were previously detected and shown in the NCBI annotation (grey
820 boxes, bottom) and are also found in our LAST alignments. NCBI annotation of genes are in green.

821 **Figure 7. Insertions and deletions in *Apis mellifera* subspecies.**
Analysis of NUMT insertions detected in only one assembly. Top: dotplot representation of LAST alignments between the two assemblies show a 745 bp variant present in AMelMel1.1 on chromosome 2 and absent in HAv3.1 (left) and a 576 bp variant present in HAv3.1 chromosome 10 and absent in AMelMel1.1 (right). For each variant, sequencing depths were evaluated on the reference in which it is present. Middle: mean sequencing depth over 80 samples (red) shows a drop coinciding with the position of the variants, suggesting that a significant proportion of samples may lack the corresponding segment and standard deviation (blue) increases in the same region, confirming the heterogeneity of the samples for the presence or absence of the variant. Bottom: mean sequencing depth per subspecies, with A. m. caucasia (15 samples) in green, A. m. ligustica (30 samples) in yellow and A. m. mellifera (35 samples) in black. Results suggest most of the A. m. mellifera samples contain the insertion present in the AMelMel1.1 assembly on chromosome 2, as the sequencing depth remains constant throughout the region, and not the one present in the HAv3.1 assembly on chromosome 10, as indicated by a sequencing depth close to zero. Inversely, most of the A. m. ligustica samples contain the insertion present in the HAv3.1 assembly on chromosome 10 and not the one in the AMelMel1.1 assembly on chromosome 2. Most A. m. caucasia samples lack the insertion present in the AMelMel1.1 assembly and a few seem to have the insertion present in the HAv3.1 assembly.

**Figure 8: Comparing the indel variant calling between sequencing depth analysis and Graphtyper2.** The Presence or absence of the NUMTs in the samples was evaluated by the pangenome graph approach with Graphtyper2 (x-axis) and by estimating the sequencing depth at the position of the NUMTs on the genome in which it is present (y-axis). Sequencing depths were normalised by calculating the ratio between sequencing depth at the position of the NUMT sequence and that of the flanking sequence. Nine out of 80 samples (11 %) could not be called for
845 NUMT_Chr2 and 19 (24 %) for NUMT_Chr10. When alleles could be called by Graphtyper2,
846 results agreed with the data based on sequencing depth.

847
Tables

<table><caption>Table 1: Literature comparison of *Apis mellifera* genetic maps</caption><thead><tr><th></th><th>Data</th><th>Physical size (Mb)</th><th>Genetic size (M)</th><th>CO/chromosome</th><th>cM/Mb</th></tr></thead><tbody><tr><td>Hunt and Page (1995)</td><td>microsatellites</td><td>178</td><td>34.5</td><td>4.3</td><td>19.4</td></tr><tr><td>Solignac et al. (2004)</td><td>microsatellites</td><td>178</td><td>40.6</td><td>-</td><td>22.8</td></tr><tr><td>Solignac et al. (2007)</td><td>microsatellites</td><td>186</td><td>40</td><td>-</td><td>22.04</td></tr><tr><td>Beye et al. (2006)</td><td>microsatellites</td><td>238</td><td>45.5</td><td>5.7</td><td>19</td></tr><tr><td>Liu et al. (2015)</td><td>SNP</td><td>220</td><td>81.4</td><td>5.1</td><td>37</td></tr><tr><td>Wallberg et al. (2015)</td><td>SNP</td><td>229</td><td>59.5</td><td>-</td><td>26</td></tr><tr><td>Wallberg et al. (2019)</td><td>SNP</td><td>219</td><td>47.3</td><td>-</td><td>21.6</td></tr><tr><td>AMelMel1.1</td><td>SNP</td><td>220</td><td>50</td><td>3.1</td><td>23</td></tr></tbody></table>
854 **Additional files**

855 **Additional file 1 Supplementary methods and Supplementary Figures S1-S18**

856 Format: pdf

857 Title: Supplementary methods and supplementary figures S1-S11

858 Description:

859 **Additional file 2 Tables S1-S11**

860 Format: Excel file

861 Title: Supplementary tables S1-S11

862 Description:

863 **Additional file 3 AMelMel-Hav3**

864 Format: pdf

865 Title: Comparison of AMelMel1.1 and HAv3.1 genome assemblies.

866 Description: Dot plot alignments of the AMelMel1.1 and HAv3.1 genome assemblies.

867 **Additional file 4 AMelMel-Amel4_5**

868 Format: pdf

869 Title: Comparison of AMelMel1.1 and Amel4.5 assemblies

870 Description: Dot plot alignments of the AMelMel1.1 and HAv3.1 genome assemblies.

871 **Additional file 5 Inversions**

872 Format: pdf

873 Title: Inversions larger than 1 kb detected between the AMelMel1.1 and HAv3.1 genome assemblies.

874

875 Description: Inversion structural variants larger than 1 kb, detected after aligning the AMelMel1.1 and HAv3.1 genome assemblies with LAST.

876
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>
# Tandem repeats
Period size: 50-1000 bp. Nb repeats: 10-1000

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>Chromosome 3</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
# 371 bp repeat

bioRxiv preprint doi: https://doi.org/10.1101/2023.12.06.570386; this version posted December 7, 2023. The copyright holder for this preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-ND 4.0 International license.

<figure><img src="image_9.png" /></figure>

# 91 bp repeat

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>
AMelMel1.1 Chromosome 3 (Mb)

<figure><img src="image_12.png" /></figure>

AMelMel1.1 Chromosome 3 (Mb)

<figure><img src="image_13.png" /></figure>

HAv3.1 Chromosome 3 (Mb)

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

Repeats identified by RepeatMasker
104, 2018-09-19
Genes, NCBI Apis mellifera Annotation Release
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>