Sonia E Eynard
email: sonia.eynard@inrae.fr

Christophe Klopp
email: christophe.klopp@inrae.fr

Kamila Canale-Tabet
email: kamila.tabet@inrae.fr

William Marande
email: william.marande@inra.fr

Céline Vandecasteele
email: celine.vandecasteele@inra.fr

Céline Roques
email: celine.lopez-roques@inrae.fr

Cécile Donnadieu
email: cecile.donnadieu@inrae.fr

Quentin Boone
email: quentin.boone@inrae.fr

Bertrand Servin
email: bertrand.servin@inrae.fr

Alain Vignal
email: alain.vignal@inrae.fr.

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

Background

The actual honey bee reference genome, HAv3.1, was produced from a commercial line sample, 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 nonprofessional apicultural practices, together with the necessity of deciphering genome backgrounds in hybrids, triggered the necessity for a specific genome assembly. Moreover, having several highquality genomes is becoming key for taking structural variations into account in pan-genome 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 a re-estimation of the honey recombination rate, over-estimated in some previous studies, due to mis-assemblies resulting in spurious inversions in the older reference genomes. The sequence continuity obtained is very high and the only limit towards continuous chromosome-wide sequences seem to be due to tandem repeat arrays usually longer than 10 kb and belonging 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, for genotyping two structural variants by a pan-genome graph approach with Graphtyper2. 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 respectively.

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 [START_REF] Sheppard | Apis mellifera pomonella, a new honey bee subspecies from Central Asia[END_REF]. The evolutionary origin of Apis mellifera is still unclear, with a possible origin in Eastern Africa or the Middle East, followed by the 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 [START_REF] Whitfield | Thrice out of Africa: ancient and recent expansions of the honey bee, Apis mellifera[END_REF][START_REF] Han | From where did the Western honeybee (Apis mellifera) originate?[END_REF][START_REF] Cridland | The Complex Demographic History and Evolutionary Origin of the Western Honey Bee, Apis Mellifera[END_REF][START_REF] Dogantzis | Thrice out of Asia and the adaptive radiation of the western honey bee[END_REF]. 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, either to be bred as pure lines or as hybrids generated by artificial or directed insemination [START_REF] Wragg | Complex population structure and haplotype patterns in the Western European honey bee from sequencing a large panel of haploid drones[END_REF][START_REF] Pieplow | A scientific note on using large mixed sperm samples in instrumental insemination of honeybee queens[END_REF]. As a consequence, these imported subspecies and hybrid lines will mate naturally to local A. m. mellifera populations, threatening them and prompting the establishment of conservation programmes [START_REF] Pinto | Genetic integrity of the Dark European honey bee ( Apis mellifera mellifera ) from protected populations: a genome-wide assessment using SNPs and mtDNA sequence data[END_REF]. However, although it has . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted December 7, 2023. ; https://doi.org/10.1101/2023.12.06.570386 doi: bioRxiv preprint been replaced in the majority of large professional beekeeper's facilities by imported honey bees, A. m. mellifera is still used by dedicated breeders.

The honey bee reference genome, whose first version was obtained as soon as 2006 [START_REF] Weinstock | Insights into social insects from the genome of the honeybee Apis mellifera[END_REF], was updated twice: a first time in 2014 [START_REF] Elsik | Finding the missing honey bee genes: lessons learned from a genome upgrade[END_REF] 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 [START_REF] Wallberg | A hybrid de novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF]. The sample used for this reference genome is from a commercial line (DH4), which is not precisely genetically defined, but is thought to be mainly of A. m. ligustica descent [START_REF] Weinstock | Insights into social insects from the genome of the honeybee Apis mellifera[END_REF]. As a consequence, the genome of the genetically distinct A. m. mellifera may not be accurately represented and future pangenome approaches, that were shown in other species to expand the number of genomic regions available for analysis [START_REF] Talenti | A cattle graph genome incorporating global breed diversity[END_REF][START_REF] Wang | The Human Pangenome Project: a global resource to map genomic diversity[END_REF], 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. This very small island (15.5 km 2 ) is located 20 km off the coast of Brittany, the conservation population was set up starting in 1987, and further imports of other honey bees banned since 1991. Mitochondrial DNA analyses have shown a low haplotype diversity and the presence of only the M-type in this population [START_REF] Garnery | Genetic diversity of the west European honey bee (Apis mellifera mellifera and A. m. iberica) I. Mitochondrial DNA[END_REF]. As expected from such a small population, microsatellite analysis has shown a low diversity [START_REF] Garnery | Genetic diversity of the west European honey bee (Apis mellifera mellifera and A. m. iberica) II. Microsatellite loci[END_REF].

Until the latest update [START_REF] Wallberg | A hybrid de novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF], the current honey bee genome sequence, Amel4.5 [START_REF] Elsik | Finding the missing honey bee genes: lessons learned from a genome upgrade[END_REF] 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 . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted December 7, 2023. ; https://doi.org/10.1101/2023.12.06.570386 doi: bioRxiv preprint 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 [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF], 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 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, 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.

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 [START_REF] Koren | Canu: scalable and accurate long-read assembly via adaptive k -mer weighting and repeat separation[END_REF] 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 [START_REF] Frith | Split-alignment of genomes finds orthologies more accurately[END_REF].

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

All the Illumina paired-end sequences from [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] 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 [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF], 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 [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF]. 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). 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 following one another within a same contig define a segment with no observed crossing over in the drones of the colony and were grouped into bins. Not having access to grand-parental 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, occurring usually on short DNA fragment often considered shorter than a few kb, [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] were removed to avoid inflating the size of the genetic map. In our study, non-crossing over gene conversion events were identified as: i) bins of length shorter than 2 kb, 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 is 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 minimising recombination events was obtained for each contig in 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. The 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 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 [START_REF] Li | Fast and accurate long-read alignment with Burrows-Wheeler transform[END_REF] using standard parameters. Variant detection was done with freebayes version 1.1.0 [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] 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 were done when alternative alleles were found in the VCF file using vcfconsensus from the vcftools package (version 0.1.12a) [START_REF] Danecek | The variant call format and VCFtools[END_REF] 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 following the same procedure as for the de-novo assembly. GC content and sequence coverage for the queens' genotypes in AMelMel1.1 were measured in 0.5Mb windows and the recombination rates were estimated using a script from [START_REF] Petit | Variation in Recombination Rate and Its Genetic Determinism in Sheep Populations[END_REF] [START_REF] Petit | Variation in Recombination Rate and Its Genetic Determinism in Sheep Populations[END_REF] over 1Mb windows. Completeness of the assemblies was estimated with BUSCO 3.0.2 [START_REF] Waterhouse | BUSCO Applications from Quality Assessments to Gene Prediction and Phylogenomics[END_REF] using OrthoDB v9.1 single-copy orthologs [START_REF] Zdobnov | OrthoDB v9.1: cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs[END_REF], 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 [START_REF] Frith | Split-alignment of genomes finds orthologies more accurately[END_REF]. Standard output psl files were produced to keep all alignments related to repeat elements, together with psl files from split alignments [START_REF] Frith | Split-alignment of genomes finds orthologies more accurately[END_REF], 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 [START_REF] Zhao | CrossMap: a versatile tool for coordinate conversion between genome assemblies[END_REF] 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) [START_REF] Benson | Tandem repeats finder: a program to analyze DNA sequences[END_REF], 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 [START_REF] Katoh | MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability[END_REF].

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 [START_REF] Waterhouse | Jalview Version 2--a multiple sequence alignment editor and analysis workbench[END_REF] 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.

Indels were detected by aligning the two genomes HAv3.1 and AMelMAl1.1 to one another with minimap2 [START_REF] Li | Minimap2: pairwise alignment for nucleotide sequences[END_REF], followed by variant calling with SVIM-asm [START_REF] Heller | SVIM-asm: structural variant detection from haploid and diploid genome assemblies[END_REF]. 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 [START_REF] Wragg | Complex population structure and haplotype patterns in the Western European honey bee from sequencing a large panel of haploid drones[END_REF] [START_REF] Wragg | Complex population structure and haplotype patterns in the Western European honey bee from sequencing a large panel of haploid drones[END_REF] and sequencing depth was estimated using SAMtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF].

Individual genotypes in the samples sequencing data was determined for the two indels by two methods. One method consisted in using GraphTyper2 [START_REF] Eggertsson | GraphTyper2 enables population-scale genotyping of structural variation using pangenome graphs[END_REF], 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 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 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 major improvement in comparison to the 46 kb N50 contig of Amel4.5 and quite similar to the N50 contig of 5.4 Mb observed in the HAv3.1 assembly [START_REF] Wallberg | A hybrid de novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF]. 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 [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF] (see Additional file 2: Table S4).

Chromosomal assignation and ordering contigs with crossing-over data.

A priori chromosomal assignment of contigs was done by alignment to the Amel4.5 assembly using LAST v956 [START_REF] Frith | Split-alignment of genomes finds orthologies more accurately[END_REF]. Out of the 200 contigs, 110 aligned successfully. Crossing-over data for confirming assignation and ordering contigs along chromosomes was obtained by using the reads from the sequencing of 43 drones from three colonies, initially used to estimate recombination rate in honey bee [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF]. Briefly, this data set contains sequence for three queens and their drone offspring (15 to 13 depending on the colony). Three of the drones of colony 1 are sequenced in duplicate and are used for the quality control of SNP calling. Aligning these reads to our contigs allowed the detection of 2,103,924 SNPs, on 176 contigs before any quality control. Out of these, approximately 64.5% were discarded due to an absence of polymorphism within the three colonies analysed, 1% for being homozygous in the queens and 1% for being heterozygous in the drones. Furthermore, 0.2% of the SNPs were discarded for being inconsistent between the three drone replicates and finally 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 reliable SNPs (approximately 32% of the initial SNPs), were detected in each of the three colonies on 114, 112 and 113 contigs respectively (see Additional file 1: Fig. 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 overall in the 104 contigs informative across all three colonies (see Additional file 1: Fig. S6). Genotype vectors for each SNP across colony drones were then defined, allowing for within-contig crossingover 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 unplaced (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 preventing assembly. To orientate the large contig on chromosome 7, positioned as 5 th 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 5 th contig, suggesting the correct orientation of the 5 th 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 [START_REF] Robertson | Canonical TTAGG-repeat telomeres and telomerase in the honey bee, Apis mellifera[END_REF] 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 were found only at the beginning of chromosome 1, whereas in the HAv3.1 assembly, they could be found at both extremities of this chromosome, but were absent from chromosomes 5 and 11 [START_REF] Wallberg | A hybrid de novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF].

An AATAT repeat was found at the beginning of chromosome 15 in our assembly.

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 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 are at the start of chromosomes 9 and 14, at the opposite end from the TTAGG repeats. On the other four chromosomes, they are at least at 1.8 Mb from a chromosome end (Fig. 2).

Recombination pattern

Having used crossing-over detection and a genetic map for contig scaffolding, we could estimate the total genetic map for AMelMel, which is approximately 50 Morgans long, giving an average recombination rate in the genome of 23 cM/Mb, close to the first estimations based on the microsatellite genetic map and to the most recent ones based on SNPs (Table 1). However, although we used the same sequencing dataset as in [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF], we found a drastic reduction in recombination rate between our genetic map and the one they initially published, which is 37 cM/Mb (Table 1). A great difference is that the latter is based on alignments of the sequence reads on Amel4.5. When aligning our assembly with 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 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 [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF]. 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 genome and 131 arrays of the 371 bp repeats. However, these numbers drop to 43 and 74 respectively when only considering tandem arrays of more than 10 repeats, suggesting that most of the 91 bp repeats have less than 10 elements (see Additional file 1: Fig. S12). To investigate sequence homogeneity within each of the two repeat classes, we selected the repeat sequence defined by TRF for repeats having strictly more than ten copies in tandem within an array. For the 91 bp repeat, we selected for 91 ≤ period size ≤ 93 and for the 371 bp repeat 367 ≤ period size 371, as suggested by the graph shown in Additional file 1: Fig. S9. Then, for each repeat class, we performed a multi-sequence alignment with MAFFT, and produced an average distance tree with Jalview. Results show 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 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 was used for a BLAST search in the AMelMel1.1 assembly. The homogeneity of the 371 bp consensus sequence was confirmed by the detection of a very high number of hits of high similarity covering the overall length of the queries (see Additional file 1: Fig. S13). On the contrary, for the three different consensus sequences used separately for the 91 bp repeat, alignment length and sequence similarities were lower, confirming that it to correspond more to a class size, rather than a specific repeat family based also on sequence composition (see Additional file 1: Fig. S13).

We then searched for the possible existence of the 371 and 91 bp repeats in other organisms.

BLAST searches with each of the 371 bp repeat consensus sequences did not allow to 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, 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 [START_REF] Vignal | Promoter sequence and chromosomal organization of the genes encoding glycophorins A, B and E[END_REF]. 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 [START_REF] Gibbons | Concerted copy number variation balances ribosomal DNA dosage in human and mouse genomes[END_REF][START_REF] Hall | Thousands of high-quality sequencing samples fail to show meaningful correlation between 5S and 45S ribosomal DNA arrays in humans[END_REF][START_REF] Ding | Genomic architecture of 5S rDNA cluster and its variations within and between species[END_REF]. 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 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 [START_REF] Cao | Chromosome-scale genome assembly of the high royal jelly-producing honeybees[END_REF]. 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 [START_REF] Castro | bric à brac (bab), a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in Drosophila melanogaster[END_REF]. 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 [START_REF] Sturtevant | The Drosophila rhomboid gene mediates the localized formation of wing veins and interacts genetically with components of the EGF-R signaling pathway[END_REF]. A more complex rearrangement involves a gene labelled as a probable nuclear hormone receptor HR38, involved in synchronizing the reproductive activity in Agrotis ipsilon [START_REF] Gassias | The insect HR38 nuclear receptor, a member of the NR4A subfamily, is a synchronizer of reproductive activity in a moth[END_REF] and in the larval-pupal transition in Leptinotarsa decemlineata [START_REF] Shen | Involvement of Leptinotarsa hormone receptor 38 in the larval-pupal transition[END_REF]. 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 [START_REF] Li | Minimap2: pairwise alignment for nucleotide sequences[END_REF] and SVIM-asm [START_REF] Heller | SVIM-asm: structural variant detection from haploid and diploid genome assemblies[END_REF]. The first one, 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 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 [START_REF] Eggertsson | GraphTyper2 enables population-scale genotyping of structural variation using pangenome graphs[END_REF], 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, for the 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 the individual genotyping. So, we first used HAv3.1 as reference and the results were a genotyping call rate of 78.7 % for NUMT_Chr2 and null for NUMT_Chr10, as the line describing its potential genotypes didn't appear 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 were also concordant with the analysis based on sequencing depth, showing that when genoyping was possible with Graphtyper2, the results were consistent. Two samples were called as heterozygotes for NUMT_Chr2, when using AMelMel1.1 as reference and were counted as "no calls", given 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, whereas successful genotyping could be obtained for samples having as little as 2X sequencing depth (see Additional file: Fig. 18). Substantially, the individual genotyping results confirm 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 have a different allele. Interestingly, NUMT_Chr2 is present in all A. m. mellifera and only two A. m. ligustica samples, and absent from all other samples, whereas NUMT_Chr10 is absent from A. m. mellifera samples and present in all but one A. m. ligustica samples and four A. m. caucasia samples (Fig. 7, Fig. 8).

Discussion

AMelMel assembly quality and comparison to other honey bee assemblies

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 [START_REF] Wragg | Complex population structure and haplotype patterns in Western Europe honey bee from sequencing a large panel of haploid drones[END_REF] 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 [START_REF] Wallberg | A hybrid de novogenome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF], 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 [START_REF] Nurk | The complete sequence of a human genome[END_REF]. 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 almost all contig boundaries investigated, long stretches of tandem repeats were found (Fig. 2). Interestingly, chromosome 16, which was obtained as a single contig, has no stretch of tandem repeats exceeding 10 kb.

Genetic maps and recombination rate in the honey bee

Having used genetic recombination data to scafold 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 [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF],

which is of the same magnitude as the latest values from [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF] and also congruent with prior values [START_REF] Beye | Exceptionally high levels of recombination across the honey bee genome[END_REF][START_REF] Solignac | A third-generation microsatellite-based linkage map of the honey bee, Apis mellifera, and its comparison with the sequence-based physical map[END_REF][START_REF] Wallberg | Extreme Recombination Frequencies Shape Genome Variation and Evolution in the Honeybee, Apis mellifera[END_REF]. 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 [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF]. On closer inspection, this higher overall recombination rate in [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF] [START_REF] Liu | Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee[END_REF], 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.

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

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 [START_REF] Kaskinova | Improved Apis mellifera reference genome based on the alternative long-read-based assemblies[END_REF].

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 [START_REF] Beye | Characterization of honeybee (Apis mellifera L.) chromosomes using repetitive DNA probes and fluorescence in situ hybridization[END_REF]. This is to date still considered as the honey bee standard karyotype structure [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF][START_REF] Robertson | Canonical TTAGG-repeat telomeres and telomerase in the honey bee, Apis mellifera[END_REF]. 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 [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF] [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF].

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 [START_REF] Robertson | Canonical TTAGG-repeat telomeres and telomerase in the honey bee, Apis mellifera[END_REF]. Taken altogether, although the current sequencing data supports the actual consensus karyotype structure, we didn't find that the AluI repeat elements [START_REF] Beye | Characterization of honeybee (Apis mellifera L.) chromosomes using repetitive DNA probes and fluorescence in situ hybridization[END_REF] 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 [START_REF] Beye | Characterization of honeybee (Apis mellifera L.) chromosomes using repetitive DNA probes and fluorescence in situ hybridization[END_REF] 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 Ava1 repeat was found at the position 5 Mb, the Ava1 elements, when detected on a chromosome, were found within 2.5 Mb of the chromosome ends, reflecting the results found on HAv3.1 [START_REF] Wallberg | A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds[END_REF]. 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 [START_REF] Wallberg | Extreme Recombination Frequencies Shape Genome Variation and Evolution in the Honeybee, Apis mellifera[END_REF] or low levels of polymorphism and recombination rates [START_REF] Wragg | Complex population structure and haplotype patterns in Western Europe honey bee from sequencing a large panel of haploid drones[END_REF]. 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.

Comparing the genomes of two honey bee subspecies. The HAv3.1 assembly is based on a sample from the DH4 line, thought to be mainly of A. m. ligustica descent [START_REF] Weinstock | Insights into social insects from the genome of the honeybee Apis mellifera[END_REF]. The comparison with our Apis mellifera mellifera AMelMel1.1 assembly allows for the detection of rearrangements occurring between these two distinct genetic types, that can't be detected through 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 crossing-over [START_REF] Hastings | Mechanisms of change in gene copy number[END_REF]. 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.

No obvious function was found for these elements to date, except for the fact that a BLAST search found that the 91 bp element shows similarity of sequence to one out of two exons of Apis dorsata and Apis florea lncRNAs, suggesting these are incomplete and consequently not active in the repeat arrays. However, the annotation of the lncRNAs in Apis dorsata and Apis florea is only based on the alignment of short reads RNA-seq. More work is needed to confirm this finding concerning the 91 bp repeat and further comparisons with other bee genomes whose sequences are underway [START_REF] Pennisi | Sequencing all life captivates biologists[END_REF] 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 possible implications in diseases [START_REF] Gibbons | Concerted copy number variation balances ribosomal DNA dosage in human and mouse genomes[END_REF][START_REF] Hall | Thousands of high-quality sequencing samples fail to show meaningful correlation between 5S and 45S ribosomal DNA arrays in humans[END_REF]. It would be interesting to see if the variations of 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. 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 [START_REF] Cao | Chromosome-scale genome assembly of the high royal jelly-producing honeybees[END_REF]. 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 [START_REF] Castro | bric à brac (bab), a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in Drosophila melanogaster[END_REF], 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 [START_REF] Sturtevant | The Drosophila rhomboid gene mediates the localized formation of wing veins and interacts genetically with components of the EGF-R signaling pathway[END_REF]. 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 [START_REF] Gassias | The insect HR38 nuclear receptor, a member of the NR4A subfamily, is a synchronizer of reproductive activity in a moth[END_REF][START_REF] Shen | Involvement of Leptinotarsa hormone receptor 38 in the larval-pupal transition[END_REF].

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 Graphtyper2. Results show that Graphtyper2 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 Graphtyper2 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 [START_REF] Sirén | Pangenomics enables genotyping of known structural variants in 5202 diverse genomes[END_REF].

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. 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, bottom) and are also found in our LAST alignments. NCBI annotation of genes are in green. LG 2_14865119_14867115

Declarations

LG 1_9031904_9033051

LG 14_5395910_5397122

LG 14_5399950_5400977

LG 1_17858200_17877100 _R_LG 1_19955057_19964380

LG 1_22345680_22350287

_R_LG 5_6159048_6163084 _R_LG 1_21563114_21564571 _R_LG 1_20842082_20853886 _R_LG 12_7163634_7170152

LG 3_4026154_4027706

LG 1_19122687_19136229

LG 1_18546954_18553935 _R_LG 11_11088037_11101425 _R_LG 11_11327297_11330583

_R_LG 12_5803634_5816532

LG 12_6253141_6262940

LG 12_6564296_6570543

_R_LG 11_12449874_12454847 _R_LG 11_12592449_12595490 _R_LG 11_12835133_12840822

_R_LG 13_7825724_7835126

LG 14_4355156_4366749

LG 2_10859245_10872752

_R_LG 8_1108674_1121039 _R_LG 11_13928907_13933762

LG 16_5947031_5954888

_R_LG 4_4658924_4672512

LG 4_4673642_4674715

LG 8_8011147_8012157

_R_LG 11_11876509_11879692

LG 11_10233084_10254059 _R_LG 16_5413238_5423304

_R_LG 11_1668221_1669318 _R_LG 3_11307002_11308339 _R_LG 7_2166894_2167997

LG 12_630917_632373

LG 3_12058679_12060307 _R_LG 5_11564585_11566224

_R_LG 7_13679455_13680597

LG 4_12458405_12459341 _R__R_LG 6_15259072_15260865 Repeats identified by RepeatMasker

FiguresFigure 1 . 7 when compared to HAv3. 1 .

 171 Figures

Figure 2 .

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

Figure 3 .

 3 Figure 3. Tandem repeats of period size 90-371 bp detected in the AMelMel1.1 assembly. The

Figure 4 .

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

Figure 5 .

 5 Figure 5. Differences in copy numbers for 5S RNA ribosomal genes. Two of the loci containing

Figure 6 . 1 .

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

Figure 7 .

 7 Figure 7. Insertions and deletions in Apis mellifera subspecies.

Figure 8 :

 8 Figure 8: Comparing the indel variant calling between sequencing depth analysis and

  . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

	bioRxiv preprint	doi:	; https://doi.org/10.1101/2023.12.06.570386	this version posted December 7, 2023.	The copyright holder for this
				8

Ethics approval and consent to participate

  International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

	bioRxiv preprint	doi:	; https://doi.org/10.1101/2023.12.06.570386	. CC-BY-ND 4.0 The copyright holder for this this version posted December 7, 2023.
	Not applicable

Table 1 : Literature comparison of Apis mellifera genetic maps 850

 1 . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

	bioRxiv preprint	doi:	; https://doi.org/10.1101/2023.12.06.570386	this version posted December 7, 2023.	The copyright holder for this

)

  . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

	371 bp repeat	_R_LG12_5803634_5816532 _R_LG4_4658924_4672512 LG1_8043072_8050375 LG1_8254379_8276064 _R_LG1_8714671_8730735 _R_LG4_6311097_6323373 LG10_4861305_4878991 _R_LG3_3577918_3588080 _R_LG15_4403367_4440377 LG10_5101150_5117639 _R_LG13_5445256_5453649 LG13_5970507_5991467 _R_LG10_7238460_7261105 LG1_17227046_17249950 _R_LG1_16991079_17005146 LG3_9243520_9263892 _R_LG13_4589794_4599483 _R_LG13_4912658_4940656 LG3_9647334_9652969 _R_LG13_3916221_3936684 _R_LG13_3614878_3648528 _R_LG13_3151098_3170763 _R_LG13_2832177_2843446 _R_LG2_5336347_5354632 _R_LG10_2029236_2044529 _R_LG13_7193095_7209886 LG13_6861162_6882002 LG1_16694845_16719536 _R_LG14_3909639_3937413 LG10_890108_919367 _R_LG10_1347023_1355959 LG10_1264396_1296486 _R_LG6_9871624_9897908 LG6_3260301_3292768 _R_LG6_2591929_2621410 _R_LG6_2495821_2505768 _R_LG6_2170060_2182944 LG2_4080811_4093590 _R_LG11_8594165_8605274 _R_LG11_8605377_8617653 _R_LG2_5981372_6004624 _R_LG2_5464450_5486877 _R_LG1_23469415_23477290 LG2_6244806_6263467 LG7_1058885_1077695 LG1_21963241_21971094 _R_LG12_6903807_6927817 LG10_1964756_1981054 LG10_1479855_1486666 LG10_1689955_1713685 _R_LG10_2526147_2541173 LG10_1173447_1184137 _R_LG10_1486767_1495772 LG10_1162734_1173346 LG10_770635_777893 LG6_7841351_7852799 LG6_7825889_7841249 _R_LG6_6217480_6240872 _R_LG6_5767962_5781061 LG6_4173483_4198584 _R_LG13_4287916_4310319 _R_LG5_9255030_9268323 LG5_3806129_3835332 LG10_1356090_1381109 _R_LG5_3671127_3688829 _R_LG6_2767784_2778126 _R_LG6_2761786_2767683 _R_LG5_3970360_3985099 _R_LG5_3985200_3992508 _R_LG5_2849370_2858886 _R_LG5_6444955_6460161 LG6_7604442_7628499 _R_LG12_8234039_8254534 LG9_8149181_8160132	91 bp repeat	Group 3 Group 2 Group 1
	bioRxiv preprint	doi:	; https://doi.org/10.1101/2023.12.06.570386	this version posted December 7, 2023.	The copyright holder for this

  . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted December 7, 2023. ; https://doi.org/10.1101/2023.12.06.570386 -seq exon coverage, aggregate (filtered), NCBI Apis mellifera Annotation Release 104 -log 2 scaled Genes, NCBI Apis mellifera Annotation Release 104, 2018-09-19 . CC-BY-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

		1,402 K	12,020 K	12,022 K	1,403 K 12,026 K 12,024 K	12,028 K	1,404 K 12,032 K 12,030 K	12,034 K	1,405 K 12,038 K 12,036 K	12,040 K	1,406 K	1,407 K	1,408 K	1,409 K	1,410 K	1,411 K	1,412 K	1,413 K	1,414 K	1,415 K	1,416 K	1,417 K	1,418 K	1,419 K	1,420 K	1,421 K	1,422 K	1,423 K
		2			2		2			2			2
		2553													2553					2553						2553
		128													128						128						128
								LOC113218699															LOC113218698	LOC113218693	LOC10...
	bioRxiv preprint	doi:						XR_003304481.1															XR_003304480.1	XR_003304475.1	exon
					LOC113218682		LOC113218684														LOC113218686	LOC113218705	LOC113218692
					XR_003304464.1	XR_003304466.1														XR_003304468.1	XR_003304486.1	XR_003304474.1
						LOC113218683			LOC113218685													LOC113218688	LOC113218689
						XR_003304465.1		XR_003304467.1													XR_003304470.1	XR_003304471.1
							LOC113218677		LOC113218695												LOC113218697	LOC113218690
							XR_003304458.1		XR_003304477.1												XR_003304479.1	XR_003304472.1
	1.380																										LOC113218691 XR_003304473.1
		12.020		12.025	12.030		12.035	12.040
					HAv3.1 Chromosome 3 (Mb)
	1.365 1.370 1.375 AMelMel1.1 chromosome 3 (Mb) this version posted December 7, 2023. ; https://doi.org/10.1101/2023.12.06.570386 doi: bioRxiv preprint	The copyright holder for this
	1.360
							1.405			1.410			1.415		1.420
														HAv3.1 chromosome 3 (Mb)

RNA

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.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests Authors' 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.

Additional files Additional file 1 Supplementary methods and Supplementary Figures S1-S18


Antoine Rouxel
email: arouxel@laas.fr

Antoine Monmayrant

Simon Lacroix

Henri Camon

Sébastien Lopez

Accurate ray-tracing optical

teaching and research institutions in France or abroad, or from public or private research centers.

The studied architecture is a single-disperser prism-based 85 CASSI (SD-CASSI), originally proposed in [START_REF] Wagadarikar | Single disperser 693 design for coded aperture snapshot spectral imaging[END_REF] and depicted in We focus throughout the paper on a prism-based SD-CASSI 96 system, but a similar approach can be trivially applied to a 97 grating-based SD-CASSI. Also, the proposed model can de-98 scribe a dual-disperser CASSI sytem (DD-CASSI), by using it 99 twice: once for each of the two arms of the DD-CASSI. 107

For the sake of simplicity, we assume the optical elements 108 and detector to be 100 % efficient: spectral filtering only occurs 109 in the mask plane. In the models proposed by [START_REF] Diaz | Adaptive multisensor 702 acquisition via spatial contextual information for compressive spec-703 tral image classification[END_REF][START_REF] Arguello | Higher-order computational 706 model for coded aperture spectral imaging[END_REF][START_REF] Gehm | Single-shot compressive 730 spectral imaging with a dual-disperser architecture[END_REF][START_REF] Dunlop-Gray | Experimental demon-733 stration of an adaptive architecture for direct spectral imaging classifi-734 cation[END_REF], H is 110 written as: 122 assumption H2: Whatever the wavelength, the magnifica-123 tion G x between the detector plane and the mask plane is 124 constant in the dispersion direction x.

125 assumption H3: Whatever the wavelength, the magnifica-126 tion G y between the detector plane and the mask plane is 127 constant in the direction perpendicular to the dispersion y.

128

The assumption H1 is contained in the function Ω(λ) that This assumption is well verified in practice when using 158 achromatic lenses and working with a small aperture system.

159

This position-to-angle relation is:

l p,a : α = -arctan (x/F) β = -arctan (y/F) ( 3
)
where F is the focal length, (x, y) the horizontal and vertical 161 incoming positions in the object plane and (α, β) the horizontal 162 and vertical outgoing angles after the lens.

163

In terms of normalized wave-vector v = k × λ 2π , after the 164 lens, the model l p,a gives:

165      v x v y v z      = 1 1 + tan 2 (α) + tan 2 (β)      tan(α) tan(β) 1      (4)
Lens model l a,p (Fourier plane to image plane) 166

Conversely, when propagating from a Fourier plane to an im-167 age plane, a perfect thin lens performs an angle-to-position con-168 version as depicted in figure 3.

Fig. 3.

Angle-to-position conversion l a,p induced by a perfect thin lens between Fourier and image planes.

169

This position-to-angle relation is:

170 l a,p : x = F tan(α) y = F tan(β) (5)
where F is the focal length, (α, β) the horizontal and vertical incoming angles and (x, y) the horizontal and vertical positions 172 in the image plane after the lens.

173

Note that using these definitions of l p,a and l a,p ensures that cation ratio of G x = G y = -1.

176

Prism model p a,a (spectral deviation)

Prisms perform an angle-to-angle conversion that depends on the wavelength. We describe the incident wavevector k i as:

k i = 2π λ × v i
where the normalized vector v i can be written as:

     v i,x v i,y v i,z      = 1 1 + tan 2 (α i ) + tan 2 (β i )      tan(α i ) tan(β i ) 1      (6)
with α and β corresponding to incident angles on the input in-

     v ox v oy v oz      =      cos(A)v ix -sin(A) [n(λ)] 2 -v 2 ix -v 2 iy v iy v o 2 -v 2 ox -v 2 oy      (7)
where A is the apex angle of the prism, v o =

208

The latter is a design choice made with respect to the spectral 209 range of the instrument.

210

The full model M p is defined by the composition of the ele-211 mentary transformations l p,a , p a,a and l a,p :

212 M p : (x d , y d , λ) → (x m , y m , λ) M p = l a,p • R P o →L 1 • p a,a • R L 2 →P i • l p,a (9)

244

An alternative is the Amici prism [START_REF] Hagen | Compound prism design principles, i[END_REF], a composite of three 245 prisms arranged to combine the benefits of both grating and 246 prism. It achieves linear dispersion over an extensive spectral 247 range, akin to gratings, while retaining the wide-range appli-248 cability of prisms [START_REF] Wu | Development of a 743 digital-micromirror-device-based multishot snapshot spectral imaging 744 system[END_REF]. Despite its merits, the Amici prism is 249 not without drawbacks, as it still manifests smile and stretch-250 ing distortions.

251

In this article, we chose to focus on the prism as it induces 252 various optical distortions that need to be addressed, and is also 253 the basic building block of the Amici prisms. Note that a grat-

254
ing model is also included in our associated simulator [START_REF] Rouxel | Simca[END_REF]. The simulated system is an SD-CASSI as depicted in figure 5.

262

We first assume that there is no misalignment (δα c = δβ c = 0) rect application of the Snell-Descartes law [START_REF] Longhurst | Geometrical and physical optics[END_REF]:

274 D m = 2 arcsin(n(λ c ) sin(A/2)) -A (10)
where n(λ c ) is the refractive index of the prism at the central 275 wavelength and A the apex angle of the prism. We finally get 276 the angle of incidence on the prism α c :

277 α c = A + D m 2 (11)
In our case, α c = 49.44°. Since the misalignments cannot be

EXPERIMENTAL VALIDATION

386

We present in this section a validation of the proposed model nique. This is more straightforward than building an experi-397 mental LUT for each individual pixel, which is reckoned as a 398 laborious process in the literature [START_REF] Wagadarikar | Video rate 723 spectral imaging using a coded aperture snapshot spectral imager[END_REF][START_REF] Poon | Practical considerations in experimental computational 726 sensing[END_REF][START_REF] Deng | Snapshot hyperspectral imaging via 748 spectral basis multiplexing in fourier domain[END_REF], even though effi-399 cient approaches have been proposed [START_REF] Ma | An efficient calibration method 751 for multi-spectral imaging[END_REF]. Also, while buidling Our SD-CASSI prototype is composed of two lenses, two 418 folding mirrors, one prism, a matrix of micro-mirrors acting as 419 a programmable mask, and a CMOS camera (figure 15). The BK7 material, equilateral, apex angle A = 60°D MD Texas Instruments DLP3010, 720 × 1280 square 5.4 µm micro-mirrors, defining a 3.9 × 6.9 mm 2 mask Detector

The Imaging Source, 33UP2000, 1920 × 1200 4.8 µm pixels

Table 1. Characteristics of the prototype components

We experimentally estimated the system point spread func- ).

An alternate way to assess the quality of the calibration is

SPECTRAL SUPER-RESOLUTION

572

In most papers [START_REF] Arguello | Higher-order computational 706 model for coded aperture spectral imaging[END_REF][START_REF] Gehm | Single-shot compressive 730 spectral imaging with a dual-disperser architecture[END_REF][START_REF] Dunlop-Gray | Experimental demon-733 stration of an adaptive architecture for direct spectral imaging classifi-734 cation[END_REF] authors assume that the data

601

The raw data-cube is then interpreted to retrieve the spectrum 602 at every point in the mask plane.

603

The figure 21 shows the spectral sampling ((a) and (c)) and

Fig. 1 .

 1 Fig. 1. Single-disperser CASSI architecture.

86 figure 1 .

 1 figure 1. Imaging optics are used to image the scene in the mask 87

  95

100A.

  Classical propagation model 101 The acquisition model widely used in the literature for an SD-102 CASSI describes the intensity collected on the detector at the 103 (x d , y d ) position as: 104 I(x d , y d ) = H(x d , y d )S 0 (x d , y d , λ) dλ (1) where H is the filtering function, (x d , y d ) the coordinates in the 105 detector plane, λ the wavelength of interest and S o the spectro-106 spatial power density of the scene.

111H 2 )

 2 (x d , y d , λ) = T(x m , y m ) with M c : x m = -x d + Ω(λ) y m = -y d (where the "classical" model M c is the backward propagation 112 model mapping positions (x d , y d ) in the detector plane to posi-113 tions (x m , y m ) in the mask plane, T(x m , y m ) is the transmission 114 efficiency of the mask at position (x m , y m ), independent of the 115 wavelength λ, and the spectral spread function Ω(λ) models 116 the spatial shift due to the spectral dispersion of the prism. 117 The underlying assumptions of the propagation model M c 118 are the following: 119 assumption H1: The spectral dispersion induced by the 120 prism is the same for all positions (x d , y d ) in the detector 121 plane.

129145

  figure 2.

Fig. 2 .

 2 Fig.2. Position-to-angle conversion l p,a induced by a perfect thin lens between object and Fourier planes.

202B. 2 .lustrated in figure 5 .

 25 M p : combination of elementary models 203 The model M p describes the propagation from a position 204 (x d , y d ) at wavelength λ through two lenses and a prism as il-205 The optical layout is such that the optical 206 axis of the lenses L 1 and L 2 corresponds to the minimum de-207 flection of the prism for the central wavelength of the system.

Fig. 5 .

 5 Fig. 5. The backward propagation model M p is the combination of the elementary transformations l p,a , p a,a and l a,p and the rotation operators R L 2 →P i and R P o →L 1 .

Fig. 6 .of the prism. 230 C. Discussions on dispersive elements 231

 6230231 Fig. 6. Prism alignment parameters. The angle α c is the angle of incidence corresponding to the minimum deflection of the prism at the central wavelength. The angles δα c and δβ c describe small errors of alignment respectively with respect to the minimum deflection and with respect to the verticality of the prism.

256

  To quantify the advantages of the model M p over the classical 257 model M c , both were compared to simulations performed by 258 the optical design software Zemax. Zemax was chosen as the 259 reference for the evaluation because it is widely used in the op-260 tical systems industry.261

Fig. 7 .

 7 Fig. 7. Typical comparison: a regular grid of 19 × 19 points in the detector plane (top) is propagated to the mask plane (bottom) with both the M c model (black dots) and Zemax (green dots).

  278 taken into account by the classical model we first assume that 279 they are non-existent. In the last validation step, we perform 280 a comparison between our model and the Zemax simulations 281 with increasing misalignments. A. Evaluation of the classical model M c 283 Preliminary methodological note: The M c model assumes iden-284 tical spectral dispersion for all positions in the field of view 285 (assumption H1, section A). The Ω function related to the disof evaluating this dispersion have been proposed in the litera-288 ture [7, 11, 14]. To avoid any bias due to the use of an inade-289 quate model for Ω, we extrapolate this function from the Zemax 290 reference simulation with a polynomial of degree 10. We obtain 291 a dispersion function for the point at the centre of the detector 292 (x d , y d ) = (0, 0). This function faithfully reproduces the spatial 293 shifts observed at the mask for the different wavelengths.

Fig. 8 .

 8 Fig. 8. Comparison between Zemax and M c propagation models for the central wavelength λ = 532 nm: (a) points in the mask plane (x m , y m ) predicted by M c (black) and by Zemax (green); (b) zoom on the stretching distortion along the horizontal axis x; (c) zoom on the smile distortion along the vertical axis y; (d) geometric distortions mapped on detector points: Euclidean distance D between points predicted by Zemax and M c .

Figure 8 (Fig. 9 .Fig. 10 .

 8910 Figure 8 (a) shows the points propagated by Zemax (black dots) and by the M c model (green dots) for the wavelength λ = 296

Fig. 11 .

 11 Fig. 11. Distance maps comparing the proposed M p model to Zemax for (left) upper, (centre) central and (right) lower wavelengths.

Figure 11 presentsFig. 12 .

 1112 Figure 11 presents the distortion maps obtained with the M p

Fig. 13 .

 13 Fig. 13. Histograms of Euclidean distances D between the Zemax simulations and the M p model, with increasing misalignment, for three wavelengths. In dotted line, the 95 th percentile of D.

387Mp

  with a prototype SD-CASSI system. The model parameters 388 are first determined with a calibration process, and the spec-389 tral sampling defined by the model is validated by observing a 390 known spectrum. 391 A. Geometric spatio-spectral calibration of an SD-CASSI sys-392 tem 393 One of the benefits of a parametric model of the image forma-394 tion is that it allows a quick calibration, using only a few im-395 ages of a known geometric pattern and an optimisation tech-396

400experimental

  LUTs is tractable for low spectral and/or spatial 401 resolution, the process becomes more difficult and time con-402 suming for high resolution imagers. 403 The proposed model-based calibration approach pertains to 404 the geometry of an SD-CASSI system, as expressed by the model 405 M p : we do not consider purely photometric phenomena (quan-406 tum efficiency and vignetting).

407Fig. 14 .

 14 Fig. 14. Calibration pattern on the DMD at the scene (a), and extracted corners in the detector plane (b).

Fig. 15 .

 15 Fig. 15. 3D model of the prototype.

FTable 2 .

 2 426 tion (PSF): it slightly changes with the position in the field of 427 view and wavelength, and its diameter is bounded by 15 µm. 428 The calibration chequerboard pattern is made up by 17 × 17 429 square regions of 30 × 30 mirrors, which defines 289 corners. 430 Three lasers are used to illuminate the chequerboard pattern, 431 with respective wavelengths 450, 532 and 780 nm. In total, 432 3 × 289 = 867 (x d , y d , x m , y m , λ) quintuplets are provided to 433 the non-linear least-square regression process to estimate the 434 parameters that define the M p model. As a reminder, these pa-435 rameters are the focal length F, the apex angle A, the angle of 436 incidence on the prism α c ensuring minimum deflection at λ c , 437 and the two angles δα c and δβ c that define the alignment errors 438 of the prism. 439 The estimated parameters are given in Table 2 with the as-440 sociated 95 % confidence interval obtained with a jackknife re-441 sampling. According to the components' specifications, the fo-442 cal lengths F and apex angle A may respectively vary by 1 % 443 and 0.1 % around their specified values of 165 mm and 60°: the 444 estimated values are within these bounds. The minimum of de-445 viation is also extremely close to the expected value. As for δα c 446 and δβ c , the estimated values are coherent with the precision of 447 the experimental alignment. In particular, as verticality of the 448 prism is easier to assess than its orientation with respect to the 449 minimum of deviation, it is not surprising that |δβ c | is smaller 450 than |δα c |. 451 [µm] A [°] α c [°] δα c [°] δβ c Calibration results: value estimated for the parameters and associated 95 % confidence interval (C.I.

452Fig.

  Fig. Experimental reprojection error D for the 3 reference wavelengths.

Fig. 17 .

 17 Fig. 17. Histogram of the reprojection error D for the 3 reference wavelengths. For each wavelength, the color dashed line corresponds to the distance of the 95 th percentile.

Fig. 18 .

 18 Fig. 18. The spectral sampling at each pixel in the detector plane depends on the back-propagation to the mask micromirror centres.

562A

  close up on the main peak at 611.5 nm of the spectrum is 563 shown in figure 20 (b): it reveals that although the measured 564 spectra are perfectly consistent and correspond to the same 565 spectral shape (give or take the acquisition noise), are 566 based on different spectral samples. 567 In this simple example of a single slit DMD pattern, recom-568 bining accurately the different lines of the detector in one single 569 spectrum with higher signal to noise ratio requires to take this 570 change of the spectral sampling into account.

Fig. 20 .

 20 Fig. 20. Spectral calibration validation across the field of view: (a) spectra measured from three different lines of the detector (y d = [310, 600, 890]) give consistent results; (b) close inspection of the spectra reveal the local spectral sampling that changes across the detector.

573

  collected by CASSI systems can be interpreted with a uni-574 form spectral sampling, thus neglecting the change observed 575 in figure 20 (b). This assumption can simplify hyperspec-576 tral cube reconstruction, as uniform spectral sampling implies 577 that acquired raw data conveniently consists of well-defined 578 monochromatic planes or images. But it leads to various issues, 579 like moiré patterns, that are particularly detrimental for high 580 resolution CASSI imagers [7]. There has been some discussion 581 regarding these issues [7-10, 12], but even the most effective 582 methods, i.e correction of the filtering function H [16], result in 583 a loss of spectral resolution by a factor of 2. Conversely, model-584 based calibration gives an easy access to the spectral local sam-585 pling, with no resolution losses and can be harnessed to achieve 586 super-resolution. 587 Indeed, considering local spectral sampling gives the abil-588 ity to achieve spectral super-resolution when several adjacent 589 pixels share the same spectrum. Such groups of pixels, dubbed 590 spectrally homogeneous spatial regions are fairly common and 591 can be harnessed to speed-up hyperspectral imaging [7, 25]. To 592 demonstrate spectral super-resolution we compare the hyper-593 spectral measurements achieved with our SD-CASSI prototype 594 with either a uniform or local spectral sampling, on a spectrally 595 homogeneous scene. 596 The acquisition scheme is a simple moving-slit sequence: we 597 successively open the columns of micro-mirrors on the DMD 598 mask and record the corresponding image on the detector. The 599 mask plane is illuminated with a large single-mode He-Ne laser 600 beam ensuring an homogeneous spectrum across many pixels.

604Fig. 21 .

 21 Fig. 21. Spectral super resolution of a He-Ne laser spectrum: the uniform spectral sampling ((a) & (b)) implies that adjacent pixels provide several intensity values for a single wavelength whereas the proposed approach ((c) & (d)) properly describes the minute changes in sampling and adjacent pixels provide spectral super-resolution. With a uniform spectral sampling (see figure 21 (a)), each of

  }, as illustrated in figure 4. Note that α i and α o belong to 180 the xOz plane and β i and β o belong to the yOz plane of the re-

	178
	179	terface of the prism P expressed in the corresponding frame
	{R P i 181
		spective frames.
		Fig. 4. Propagation of a wavevector through a prism. k i is an
		incident wavevector on the input interface of the prism in the
		{R P i } frame; k o is the corresponding output vector in frame
		{R P o }.
	182
	183	The refraction law ensures that the component of the
	184	wavevector lying in the air/glass interface (i.e. component
	185	within the z = 0 plane in {R P i }) is conserved. The same applies
	186	to its normalized vector v i . Combining this conservation of in-
	187	plane components at both input and output interfaces with the
	188	fact that the normalized wavevector norm is 1 out of the prism
	189	and n(λ) inside the prism, one can express the components of
	190	the output vector v o in the {R P o } frame as:

a 4F telescope, described by l p,a • l a,p has the expected magnifi-

To evaluate the spectral sampling of an SD-CASSI or a DD-476 CASSI system, one usually resort to a an experimental LUT-477 based calibration [START_REF] Arguello | Higher-order computational 706 model for coded aperture spectral imaging[END_REF][START_REF] Wagadarikar | Video rate 723 spectral imaging using a coded aperture snapshot spectral imager[END_REF][START_REF] Poon | Practical considerations in experimental computational 726 sensing[END_REF][START_REF] Arce | Compressive coded aper-753 ture spectral imaging: An introduction[END_REF][START_REF] Goldstein | DMD-based adap-756 tive spectral imagers for hyperspectral imagery and direct detection[END_REF]. The proposed method relies