Héloïse Bastide
email: heloise.bastide@universite-paris-saclay.fr

Hélène Legout

Noé Dogbo

David Ogereau

Carolina Prediger

Julie Carcaud

Jonathan Filée

Lionel Garnery

Frédéric Clément Gilbert

The genome of the blind bee louse fly reveals deep convergences with its social host and illuminates Drosophila origins

Keywords: parasitism, inquilinism, phylogenomics, horizontal transposon transfer, gene family evolution, adaptation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Results and Discussion

The bee louse fly Braula coeca is an aberrant member of the Drosophilidae Among the several parasites and inquilines that are attracted by the rich resources and clean and protective shelter of the Western honey bee Apis mellifera nest, none has undergone as profound morphological changes as the apterous and quasi-blind bee louse fly Braula coeca (Figure 1A-C). The female lays eggs in honey (not brood) cells, and the hatched larvae eat pollen and wax, where they burrow tunnels in which they pupate without forming true puparia. [START_REF] Skaife | On Braula Coeca , Nitzsch, a Dipterous parasite of the honey bee[END_REF][START_REF] Imms | On Braula coeca Nitsch and its affinities[END_REF] Following emergence, the adults attach to the body of worker bees, migrating from one individual to another until reaching the queen (Figure 1A). There, they move to the queen's head, stimulate regurgitation, and imbibe from her mouth honey and nectar. [START_REF] Skaife | On Braula Coeca , Nitzsch, a Dipterous parasite of the honey bee[END_REF][START_REF] Imms | On Braula coeca Nitsch and its affinities[END_REF] The bee louse fly is considered an inquiline kleptoparasite with potential negative effects on honey bee colony health due to the galleries it makes in bee combs and the facilitation of transmitting serious pathogenic viruses to the bees. [START_REF] Avalos | Honey bee louse (Braula schmitzi) as a honey bee virus vector[END_REF] Ever since Réaumur's first description of the bee louse fly in 1740, and Nitzsch's creation of the genus Braula in 1818, [START_REF] De Réaumur | Mémoires pour servir à l'histoire des insectes[END_REF][START_REF] Nitzsch | Die Familien und Gattungen der Thierinsekten (insecta epizoica); als Prodromus einer Naturgeschichte derselben[END_REF] the positioning within the Diptera of the bee louse fly and affiliated species that were classified under the family Braulidae has been puzzling because of its aberrant morphology and unique adaptations to a social host. This family contains seven species belonging to the genera Braula and Megabraula that are all inquilines to honey bee species of the genus Apis. Recent phylogenetic analyses based on a transcriptome assembled from one adult fly and using 1,130 loci interestingly showed Braula coeca, the most widespread braulid, to constitute a basal lineage within the Drosophilidae that was sister to four genera of the subfamily Steganinae. [START_REF] Bayless | Beyond Drosophila: resolving the rapid radiation of schizophoran flies with phylotranscriptomics[END_REF][START_REF] Winkler | Phylogenetic resolution of the fly superfamily Ephydroidea-Molecular systematics of the enigmatic and diverse relatives of Drosophilidae[END_REF] To reassess this hypothesis using a larger dataset, we sequenced the whole genome from a pooled sample of 15 unsexed B. coeca flies, all collected on Ouessant Island in Western France.

We used a hybrid approach to assemble a genome using long-read Oxford Nanopore Technology (ONT) and short-read Illumina sequencing (see Methods). Benchmarking Universal Single-Copy Orthologs (BUSCO) [START_REF] Simão | BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[END_REF] gave a score of 95.8% of the Dipteran conserved single-copy orthologs with 1.3% of duplicated genes. This value is higher than the recommended score of 90% for reference genomes [START_REF] Lewin | Earth BioGenome Project: Sequencing life for the future of life[END_REF] . Merqury 15 estimated an assembly completeness of 93.6% and a consensus quality value (QV) of 41, which exceeds the recommended threshold of QV40 for reference genome. [START_REF] Koren | Reply to 'Errors in long-read assemblies can critically affect protein prediction[END_REF] We assembled two genomes and one transcriptome of three additional steganine genera. We then built a supermatrix of 3,100 BUSCO genes (2,557,349 amino acids) that included 15 drosophilid species (representative members of the four main radiations in the family), [START_REF] Yassin | Phylogenetic classification of the Drosophilidae Rondani (Diptera): the role of morphology in the postgenomic era[END_REF] and 5 species belonging to the superfamily Ephydroidea to which both the Drosophilidae and Braulidae belong (Table S1). The maximum-likelihood phylogenetic analysis of this large dataset reconfirmed the close-relatedness of B. coeca to the Drosophilidae. It further showed that it is a full member of the subfamily Steganinae (Figure 1D). The taxonomic priority principle should consider the family Drosophilidae, described in 1856, [START_REF] Rondani | Dipterologiae Italicae prodromus: Genera Italica Ordinis Dipterorum (Stocchi)[END_REF] a junior synonym for the family Braulidae, described in 1853. [START_REF] Egger | Himmlische Waffenrüstung für die Jugend bestehend aus den heiligen Sakramenten der Busse, des Altars und der Firmung: Ein praktischer Unterricht[END_REF] However, the asymmetric size and scientific relevance of the two families argue against such a decision. We, therefore, opt for synonymizing the Braulidae with the Drosophilidae, referring hereafter to Braula and Megabraula as members of the subfamily Steganinae.

Inquilinism in the bee louse fly likely evolved from sap breeders associated with scale insects

To gain further insight into the history of the association between Braula and Apis, we mapped the predominant ecological habitats of ephydroid families on the phylogeny.

The ancestral habitat of ephydroids was presumed to be rotting leaf molds. [START_REF] Throckmorton | The phylogeny, ecology, and geography of Drosophila[END_REF] From this, multiple specializations took place, including the exploitation of aquatic molds (and eventually algae) in the Ephydridae, [START_REF] Keiper | Biology and ecology of higher Diptera from freshwater wetlands[END_REF] Mammal dung in Curtonotidae and Diastatidae, [START_REF] Keiper | Biology and ecology of higher Diptera from freshwater wetlands[END_REF][START_REF] Pollock | Observations on the biology and anatomy of Curtonotidae (Diptera: Schizophora)[END_REF] and fermenting vegetables and fruits, sap and fungi, with specialization mostly on yeasts in the Drosophilidae 23 (Figure 1D). Bayesian reconstruction suggests the ancestral habitat of the Drosophilidae to be tree sap breeding (Figure 1D) with fungus-and fruit-breeding subsequently deriving and predominating in the genera Leucophenga and Drosophila, respectively. Remarkably, the deepest branches in the Steganinae and the Cryptochaetidae (the closest relative to the Drosophilidae) represent lineages whose larvae are predatory of scale insects and mealy bugs, e.g., Acletoxenus formosus and A. indicus on aleyrodoids, Rhinoleucophenga brasiliensis and R. obesa as well as Cryptochaetidum iceryae and C. grandicorne on coccoids. [START_REF] Ashburner | Entomophagous and other bizarre Drosophilidae[END_REF] In those lineages, adults are often seen to feed on the honeydew produced by the bugs, an abundant sugar-rich substrate sucked from plants' sap, while larvae take shelter and develop in the waxy secretions of these insects. This dependence on sugary substrate (honeydew) and development in a waxy environment could have predisposed Braula's inquilinism in bee nests.

The bee louse fly inquilinism is relatively recent

To date Braula inquilinism, we inferred a fossil-calibrated phylogeny using 79 single-copy orthologs (63,192 amino acids) in 17 Acalyptrate dipteran and 25 Apocrite hymenopteran species (see Methods; Table S1; Figure 2A). Five non-ephydroid dipteran species with Ref-Seq assemblies were included in this analysis to correct for tree imbalance. [START_REF] Duchêne | Tree imbalance causes a bias in phylogenetic estimation of evolutionary timescales using heterochronous sequences[END_REF] The divergence between B. coeca from its closest steganine relatives (node 1 in Figure 2A: 44.9 [37.8-53.8] million years (myr) [95% confidence interval]) overlapped with the origin of the Apidae (50.12 [42.9-64.5] myr) and with the transition from solitary to subsocial (2: 40.25 [32.4-47.8] myr) and primitively social habits (3: 30.98 [26.3-36.0] myr). [START_REF] Shell | Sociality sculpts similar patterns of molecular evolution in two independently evolved lineages of eusocial bees[END_REF] It is possible that the origin of the bee louse fly-apid interactions occurred at sap breeding sites, when early subsocial apids started to gather resin and other plant exudates, as well as scale insects honeydew, and stored them in their nests. As eusociality evolved (4: 23.8 [18.9-28.2] myr), the proportion of resin to secreted wax diminished, and some cells were also used to store nectar and honey for the brood. [START_REF] Noll | Behavioral phylogeny of corbiculate Apidae (Hymenoptera; Apinae), with special reference to social behavior[END_REF] A shift from the putatively ancestral dependence on honey and wax produced by scale insects to those produced by bees might have evolved by then. The transition to eusociality in the genus Apis required an important division of labor that involved the evolution of pheromonal control of the reproductive capacity of worker females by the queen and the evolution of trophallaxis. [START_REF] Noll | Behavioral phylogeny of corbiculate Apidae (Hymenoptera; Apinae), with special reference to social behavior[END_REF] Adaptation of Braula to the queen pheromone compounds that have anti-ovarian effects on a wide range of insects, including Drosophila melanogaster [START_REF] Galang | Analysis of the Drosophila melanogaster anti-ovarian response to honey bee queen mandibular pheromone[END_REF] and the exploitation of trophallaxis 7 could not have evolved before the advancement of eusociality (5: 17.1 [11.6-23.0] myr). The evolution of blindness and apterism should have constrained the dispersal of the bee lice, relating their speciation history to that of their hosts. Indeed, only seven bee louse fly species are known, of which five Braula species are restricted to the Western honey bee A. mellifera, and two Megabraula species are restricted to the giant honey bee A. laboriosa in the Himalayas. [START_REF] Grimaldi | Megabraula, a new genus for two new species of Braulidae (Diptera), and a discussion of braulid evolution[END_REF][START_REF] Dobson | A "bee-louse" Braula schmitzi örösi-pál (Diptera: Braulidae) new to the British Isles, and the status of Braula spp[END_REF] The divergence between these Apis species, and presumably between Braula and Megabraula, is estimated at 6: 5.8 [2.8 -12.0] myr ago. Therefore, the evolution of the bee louse fly inquilinism has likely taken place during the Mid-to Late Miocene period between 5.8 and 17.1 myr ago (Figure 2A). We cannot rule out even a more recent origin if the ancestor of Braula or Megabraula has shifted from one Apis host to another, i.e. <5.8 myr ago.

The bee louse fly inquilinism was accompanied by a reduction of gene content but not genome size

Loss of significant portions of genomic and gene contents is a characteristic of obligate parasites specializing on specific hosts or inhabiting extreme environments. For example, the human body louse, Pediculus humanus, has one of the smallest genomes and the lowest numbers of genes in insects (108 megabases [Mb] and 10,773 protein-coding genes). [START_REF] Kelley | Compact genome of the Antarctic midge is likely an adaptation to an extreme environment[END_REF] For B. coeca, we obtained a final assembly size of 309.35 Mb shared by 2,477 contigs with an N50 of 347,211 bp. This N50 estimate is typical of hybrid genome assemblies obtained using a pooled sample of wild-caught drosophilid flies from species with large genome sizes (> 300 Mb). [START_REF] Kim | Highly contiguous assemblies of 101 drosophilid genomes[END_REF] No evidence for polyploidy or other endosymbiont that could have biased the genome size estimate was detected (Figure S1). Genome size prediction using k-mers distribution spectra predicted a genome of 308 Mb, concordant with the assembly size (Figure S1). Such a genome size is significantly larger than the remaining drosophilid species (Student's t one-sample test P < 1. S1).

To determine the number of protein-coding genes, we used four rounds of Maker [START_REF] Cantarel | MAKER: An easy-to-use annotation pipeline designed for emerging model organism genomes[END_REF] supported by the training of the gene finding and prediction tools SNAP [START_REF] Korf | Gene finding in novel genomes[END_REF] and Augustus [START_REF] König | Simultaneous gene finding in multiple genomes[END_REF] .

The annotation, made on the repeat-masked genome, yielded 10,349 protein-coding genes with an Annotation Edit Distance (AED) ≤ 0.5 for 96.4% of our gene models and a Pfam domain found in 83.66% of the proteins (BUSCO score = 91%). Using the same strategy, we annotated two steganine genomes, namely Phortica variegata and Leucophenga varia. The annotation yielded 11,067 (BUSCO score = 91%) and 13,160 (BUSCO score = 90.8%) protein-coding genes, respectively. The annotation of the ephydrid Ephydra gracilis genome yielded 9,154 protein-coding genes (BUSCO score = 68.9%) (Figure S2). Ephydra is particular among Ephydroidea in adapting to hypersaline waters and associated algal flora. [START_REF] Van Breugel | Superhydrophobic diving flies (Ephydra hians) and the hypersaline waters of Mono Lake[END_REF] Given the current low knowledge of ephydrid genetics, whether their low gene content is due to their high specialization or an artifact of incomplete annotation, is hard to know. Regardless, the bee louse fly has the lowest number of protein-coding-genes compared to other drosophilids (Student's t one-sample test P < 1.5 x 10 -5 ; Shapiro-Wilk normality test P = 0.48) despite having a total genome size that is among the largest genomes in the family. Whereas the near-completeness of our B. coeca genome (~95%) might have reduced the number of annotated genes, low complexity hard-to-assemble genomic regions are usually mostly heterochromatic and poor in genes, e.g., centromeres, Y chromosomes, etc. A low gene content is also characteristic of bee genomes, compared to ants and wasps, with a remarkable trend of gene reduction within the family Apidae during the evolution of the genus Apis (Figure 2C; Figure S2; Table S1).

Transposable elements (TEs) expanded in the bee louse fly with one element horizontally transferred with the host

The bee louse fly's large genome size and low gene content suggest an increase in repetitive sequences. RepeatModeler and RepeatMasker analyses [START_REF] Flynn | RepeatModeler2 for automated genomic discovery of transposable element families[END_REF]38 indicated that nearly 41.34% of the B. coeca genome consists of such sequences, compared to 22.05% and 10.98% in D. melanogaster and A. mellifera, respectively (Figure 2D). Remarkably, half of the bee louse fly repetitive sequences consisted of long interspersed nuclear elements (LINEs) retrotransposons (14.94%). While LINEs are usually among the most abundant transposable elements after LTRs within the Drosophilidae, 39 their values did not exceed what was found in B. coeca (we found the highest percentage in Leucophenga varia with 5.54%). It is at present unclear what factors influence the diversity of transposable elements (TEs) landscapes among eukaryote species [START_REF] Gilbert | Transposable elements and the evolution of insects[END_REF] . Nonetheless, this difference means that whereas the bee louse fly has likely retained the ancestral large genome size of the Drosophilidae, its TEs constitution has largely evolved.

Because host-parasite relationships have repeatedly been invoked as a factor that may favor horizontal transfer of TEs, [START_REF] Venner | Ecological networks to unravel the routes to horizontal transposon transfers[END_REF] we searched for evidence of such transfers between B. coeca and A. mellifera. We found one TE, a DNA transposon Famar1-like element previously described in the earwig Forficula auricularia [START_REF] Barry | A Bacterial Genetic Screen Identifies Functional Coding Sequences of the Insect mariner Transposable Element Famar1 Amplified From the Genome of the Earwig, Forficula auricularia[END_REF] that belongs to the Tc1mariner superfamily. This element showed a high similarity between B. coeca and A. mellifera but was absent in all other drosophilid species for which a genome is available in GenBank, which is highly suggestive of an acquisition through horizontal transfer (Figure S2). Indeed, phylogenetic analysis of multiple copies of this TE extracted from 37 widely divergent animal species (Figure 2E) supported a direct transfer event between B. coeca and A. mellifera, although the directionality of the transfer cannot be inferred since the elements from the two species form mutually-exclusive monophyletic clades.

Remarkably, all elements found in the genomes of four A. mellifera subspecies, including A. m. carnica, A. m. caucasia, A. m. mellifera and A. m. ligustica, formed an exclusively monophyletic clade. The transfer time between B. coeca and A. mellifera has likely preceded the dispersion of this element among the subspecies or even their differentiation 0.77 myr ago [START_REF] Carr | Multiple mitogenomes indicate things fall apart with Out of Africa or Asia hypotheses for the phylogeographic evolution of honey bees (Apis mellifera)[END_REF] if the element was ancestral in A. mellifera. On the other hand, we did not find any trace of this element or any other related element in any other Apis species, indicating that the maximal time of horizontal transfer likely does not surpass 2.73 [0.70-8.9] myr ago, i.e. the time of divergence between A. mellifera and its closest-relative A. cerana (Figure 2A). The tight ecological connection between the bee louse fly and its host may have favored this transfer, as was suggested for blood-or sapsucking insects. [START_REF] Gilbert | A role for host-parasite interactions in the horizontal transfer of transposons across phyla[END_REF][START_REF] Gilbert | Multiple Horizontal Acquisitions of Plant Genes in the Whitefly Bemisia tabaci[END_REF]

Gene families with excess losses show striking cross-order parallelism

Despite their deep divergence, we tested whether parallel changes could explain the reduction of protein-coding genes in both the honey bee and the bee louse fly. We used OrthoFinder [START_REF] Emms | OrthoFinder: phylogenetic orthology inference for comparative genomics[END_REF] to cluster orthologous proteins from the 25 hymenopteran and 17 dipteran species. We identified 19,010 orthogroups. Of these, 935 showed significant size evolution among the 42 species when analyzed using CAFE5 [START_REF] Mendes | CAFE 5 models variation in evolutionary rates among gene families[END_REF] and after applying an error model that accounted for misassemblies and misannotations. To classify those orthogroups into functional categories, we extracted groups that contained D. melanogaster orthologs for which a molecular function, i.e. a gene group, was assigned in the Flybase database [START_REF] Thurmond | FlyBase 2.0: the next generation[END_REF] (see Methods). Of 1,078 gene groups, 136 significantly deviated from the birth-death model estimated by CAFE5.

After correction for multiple testing, 17 gene groups had significant losses in the bee louse fly with no group showing significant gain (Table 1). The reduction of most of these groups showed a striking parallelism with bees (Anthophila) in particular and hymenopterans in general (Table 1, Figure S3). The most significant groups were those involved in the chemical detection of taste (gustatory receptors and divergent ionotropic receptors) and odors (odorant receptors and odorant binding proteins). The remaining groups included those involved in recognition and signaling with a potential role in metabolism, immunity, and/or development such as C-type lectins, serine proteases, and Dorsal, [START_REF] Emms | OrthoFinder: phylogenetic orthology inference for comparative genomics[END_REF] as well as ion and sugar transportations. Other groups are involved in detoxification, such as cytochrome P450, GST-C, and carboxylases. [START_REF] Mendes | CAFE 5 models variation in evolutionary rates among gene families[END_REF] Indeed, bees have evolved a reduced repertoire of immunity and detoxification genes, likely due to the evolution of social behavior and their life in an overprotective and clean shelter, i.e. the nest. [START_REF] Evans | Immune pathways and defence mechanisms in honey bees Apis mellifera[END_REF][START_REF] Berenbaum | Xenobiotic detoxification pathways in honey bees[END_REF] Cytochrome P450 genes are more expressed in foraging workers than in the castes that remain in the nest (i.e. the queen and nurse workers). [START_REF] Chan | Honey bee protein atlas at organ-level resolution[END_REF] The reduction of peptidases in both the honey bee and the bee louse fly could also be due to the low protein content of some of their food, i.e. nectar and honey. We also noted an underrepresentation of chitin-binding domain proteins and chitinases in the bee louse fly and the honey bee.

Cuticles could act as barriers against environmental toxins, which may not be highly encountered in the nest. Remarkably, B. coeca is unique among cyclorrhaphan dipterans as its pupa, similarly to the honey bee's, [START_REF] Winston | The Biology of the Honey Bee[END_REF] is contained in the unmodified cuticle of the third instar larva, and no sclerotized puparium is formed. [START_REF] Skaife | On Braula Coeca , Nitzsch, a Dipterous parasite of the honey bee[END_REF][START_REF] Imms | On Braula coeca Nitsch and its affinities[END_REF] Whereas assembly and annotation errors can bias general estimates of gene losses, they should not specifically target the gene families that are ecologically relevant to both the host and the inquiline.

Honey and wax feeding drove the loss of almost all bitter-tasting gustatory receptors

The two most significantly evolving gene families in the bee louse fly, i.e. gustatory receptors (GRs) and divergent ionotropic receptors (IR-DIVs), allow the detection of soluble cues (Table 1). There are 60 GRs in D. melanogaster, of which 9 and 49 receptors respond primarily to sweet and bitter tastes, respectively, and 2 receptors respond to carbon dioxide (CO2). [START_REF] Weiss | The Molecular and Cellular Basis of Bitter Taste in Drosophila[END_REF] The three categories clustered into 35 orthogroups (Figure 3A), whose phylogenetic analysis indicates that the ancestral drosophilid repertoire consisted of 6 sweet, 26 bitter, and 2 CO2 GRs assuming functional conservation of gustatory categories (Figure 3A). We identified 11 GRs in the bee louse fly with no duplications using InsectOR 54 and manual curation. These GRs could be classified according to their D. melanogaster orthologs into 3 sweet, 6 bitter, and 2 CO2. That means that the D. melanogaster lineage disproportionally evolved more bitter receptors from the ancestral repertoire, whereas B. coeca disproportionally lost bitter receptors (Figure 3A). InsectOR inferred the number of GRs in the steganine species L. varia and P. variegata to be 21 and 26, respectively, further confirming that B. coeca has lost a significant portion of the ancestral GR repertoire (Figure S4). Honey bees have only 11 GRs, of which 7 are orthologous to sweet Drosophila GRs. [START_REF] Sadd | The genomes of two key bumblebee species with primitive eusocial organization[END_REF] This is likely due to the bees' strong diet reliance on sweet floral nectars and honey. [START_REF] Robertson | The chemoreceptor superfamily in the honey bee, Apis mellifera: expansion of the odorant, but not gustatory, receptor family[END_REF] The loss of B. coeca bitter GRs and its retention of 2 ancestral sweet receptors is a strong convergence with its host. Ionotropic receptors are another major class of chemoreceptors. They are divided into antennal IRs, which are conserved across insects and are most likely involved in olfaction, and divergent IRs (IR-DIVs), which evolve rapidly and are mostly involved in the taste perception of carboxylic and amino acids. Only divergent IRs showed a significant loss in B. coeca (Table 1). However, our knowledge about the function of the 42 D. melanogaster IR-DIVs is still limited. [START_REF] Ni | The structure and function of ionotropic receptors in Drosophila[END_REF] We inferred the ancestral IR-DIV drosophilid repertoire to contain 29 receptors, of which only 9 were retained in B. coeca. Remarkably, whereas we found almost no direct orthologs between Diptera and Hymenoptera for IR-DIVs (Figure S3), bees are known to have few IRs in general [START_REF] Park | Uncovering the novel characteristics of Asian honey bee, Apis cerana, by whole genome sequencing[END_REF] pointing to another possible taste convergence between the bee louse fly and its host.

One-fifth of ancestral odorant receptors was lost, including one receptor that is involved in anti-ovarian response in Drosophila melanogaster

Odorant receptors (ORs) are essential to detect volatile chemical cues from the environment. This family has expanded in the honey bee to reach 170. [START_REF] Karpe | Identification of Complete Repertoire of Apis florea Odorant Receptors Reveals Complex Orthologous Relationships with Apis mellifera[END_REF] However, only 9 of the honey bee genes have orthologs with D. melanogaster, and phylogenetic analysis indicates that this common OR repertoire has been gradually reduced during the evolution of Apis (Table 1, Figure S3). The 60 ORs of D. melanogaster clustered within 16 orthogroups (Figure 3B). We inferred the ancestral drosophilid OR repertoire to contain 44 ORs, with at least one representative for each orthogroup (Figure 3B). We identified in B. coeca, following InsectOR 54 and manual curation, 35 ORs in addition to Orco, i.e. onefifth of the ancestral repertoire was lost. The number of ORs was 50 and 51 in the two closely related steganine species L. varia and P. varia, respectively (Figure S4). Braula ORs were direct orthologs to 18 genes in D. melanogaster (Figure 3B). Judging from the response of these orthologs to different volatiles in D. melanogaster as curated in the DoOR database [START_REF] Münch | DoOR 2.0 -Comprehensive mapping of Drosophila melanogaster odorant responses[END_REF] and assuming potential conservation of function, the retained bee louse fly ORs may respond to compounds produced by honey bee workers in a defense context (e.g., 1-hexanol, farnesol, 2-heptanone), [START_REF] Wang | Honey Bee Alarm Pheromone Mediates Communication in Plant-Pollinator-Predator Interactions[END_REF] and/or of floral, pollen and nectar aromas, such as acetophenone and benzaldehyde, a major volatile of honey. [START_REF] Machado | Honey Volatiles as a Fingerprint for Botanical Origin-A Review on their Occurrence on Monofloral Honeys[END_REF][START_REF] Starowicz | Characterizing the Volatile and Sensory Profiles, and Sugar Content of Beeswax, Beebread, Bee Pollen, and Honey[END_REF] Two cases of tetraplications were observed. One case involved three recent duplications of genes orthologous to DmOr67b, a gene that is highly responsive in D. melanogaster to both acetophenone and 1-hexanol. The second case involved three successive duplications of a gene orthologous to DmOr74a, which responds in D. melanogaster larvae to 1-nonanol and 1-heptanol, the latter being a major brood volatile, 64 and 1-hexanol, a component of the alarm pheromone. [START_REF] Collins | Alarm responses caused by newly identified compounds derived from the honeybee sting[END_REF] Of these three duplications, two were unique to B. coeca compared to its closely-related steganine species (Figure S4). Low concentrations of isopentyl acetate, the main component of the alarm pheromone, released by unstressed workers at the nest entrances attract the parasitic nest beetle Aethina tumida, [START_REF] Torto | Multitrophic interaction facilitates parasite-host relationship between an invasive beetle and the honey bee[END_REF] suggesting that the detection of the host odors could be a common strategy among phylogenetically distant inquilines and parasites of social insects.

Whereas major molecular convergences could exist between the inquiline and its social host, divergent strategies to adapt to the eusocial lifestyle requirements are still needed. In honey bees, colony cohesion is driven by the volatile queen's mandibular pheromone (QMP), which "sterilizes" the bee workers. [START_REF] Van Oystaeyen | Conserved Class of Queen Pheromones Stops Social Insect Workers from Reproducing[END_REF] This pheromone elicits an antiovarian response in other insects, including D. melanogaster. [START_REF] Galang | Analysis of the Drosophila melanogaster anti-ovarian response to honey bee queen mandibular pheromone[END_REF] An RNA interference (RNAi)-screen identified DmOr49b, DmOr56a, and DmOr98a to be potentially involved in the detection of the QMP compounds and the suppression of fecundity. [START_REF] Galang | Analysis of the Drosophila melanogaster anti-ovarian response to honey bee queen mandibular pheromone[END_REF][START_REF] Camiletti | A novel screen for genes associated with pheromone-induced sterility[END_REF] A sine qua non condition for a drosophilid to reproduce in a bee nest would, therefore, be to lose those receptors or to modify their response or effect. We found that the bee louse fly does not have an ortholog for DmOr98a, a receptor specific to the genus Drosophila (Figure S4).

The bee louse fly has a pseudogene, orthologous to DmOr49b, that InsectOR identified.

Orthologs of this D. melanogaster receptor are present and complete in all dipteran species, including L. varia and P. variegata (Figure S4). The bee louse fly had a receptor that we called BcOr22, which was orthologous to DmOr56a (Figure 4B). This last receptor is narrowly tuned in many Drosophila species to a single component, the mold volatile geosmin, whose perception also inhibits oviposition in D. melanogaster, [START_REF] Stensmyr | A Conserved Dedicated Olfactory Circuit for Detecting Harmful Microbes in Drosophila[END_REF]

Blindness and life in a dark nest were accompanied by the loss of multiple rhodopsins

The species Latin name of the bee louse fly refers to the assumption that it was blind due to the reduction of the eye size and the loss of the ocelli. In agreement with reduced vision in the bee louse fly, we found only two out of the seven rhodopsin genes, which are responsible for colored vision and positive phototaxis in D. melanogaster and which were all present in the ancestral drosophilid repertoire. D. melanogaster orthologs of the Rh1 and Rh6 genes are expressed in the ommatidia and are sensitive to light. [START_REF] Senthilan | Rhodopsin 7-The unusual Rhodopsin in Drosophila[END_REF] The role of these opsins in light detection, despite the absence of ommatidia in the bee louse fly is unclear. Remarkably, Rh1 and Rh6 are structurally required in mechanosensory bristles to control larval locomotion. [START_REF] Zanini | Proprioceptive Opsin Functions in Drosophila Larval Locomotion[END_REF] They also detect temperature. [START_REF] Leung | Unconventional Roles of Opsins[END_REF] Therefore, the retention of these rhodopsins in the bee louse fly could mainly be due to their unconventional functions. On the other hand, the rhodopsin Rh2, which is exclusively expressed in the ocelli and used for horizon detection in D. melanogaster [START_REF] Mishra | Homothorax controls a binary Rhodopsin switch in Drosophila ocelli[END_REF] , is among those lost in the bee louse fly, in agreement with the loss of the ocelli. Regression of the visual system and its underlying opsin genes is common in animals inhabiting dark environments, such as fossorial mammals [START_REF] Partha | Subterranean mammals show convergent regression in ocular genes and enhancers, along with adaptation to tunneling[END_REF] and cavefishes, 75 representing a major example of deep convergences.

Apterism was not accompanied by the loss of major wing development genes

Small size, loss of wings, and the evolution of strongly clinging legs are all morphological changes that could prevent the honey bees from getting rid of the bee lice. [START_REF] Büscher | The exceptional attachment ability of the ectoparasitic bee louse Braula coeca (Diptera, Braulidae) on the honeybee[END_REF] All these potential adaptations are convergent with ectoparasitic true lice, and for some, such as apterism, represent major recurrent changes that have responded to distinct pressures throughout the history of insects. [START_REF] Roff | The Evolution of Flightlessness in Insects[END_REF] We found intact most of the main wing development genes whose mutations severely reduce the wing in D. melanogaster, such as wingless, apterous, or vestigial. This means that the major morphological changes more likely resulted from regulatory changes of these core genes or modifications of other genes. Future developmental studies, specifically comparing the expression of wing and leg morphogenic genes between the bee louse fly and D. melanogaster, will definitively help shed light on the transcriptomic shifts underlying the major morphological changes of the bee louse fly.

Conclusion

That the enigmatic bee louse fly is indeed a drosophilid, a lineage within the most investigated insect family with more than 150 fully sequenced genomes, is undoubtedly one of the most exciting discoveries in dipteran phylogeny. How could a fly with an ancestral drosophilid genome become ecologically adapted to bees and morphologically similar to lice? Our results show that a mosaic of deep convergences at the genomic level underlies the relatively recent and dramatic changes of the bee louse fly to nest inquilinism. This mosaicism involved deep convergences with the host, mostly in genes likely involved in immunity, detoxification, and chemical perception, as well as convergences with general features of fossorial animals in the visual systems. Future developmental studies may elucidate whether general morphologies, such as apterism and leg modifications, could also be shared between Braula and other ectoparasites. Due to its genetic relatedness to Drosophila and ecological association to Apis, two major laboratory models, the new genomic resources presented here can help establish the bee louse fly as a promising model to address questions related to deep convergences that are difficult to approach in multiple highly specializing animals.

Figures titles and legends

STAR Methods

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Héloïse Bastide (heloise.bastide@universite-parissaclay.fr).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• Raw sequence data are deposited on NCBI Sequence Read Archive (SRA).

Bioproject accession number is listed in the key resources table. Genome assemblies and all data associated to this study including translation of early taxonomic literature are deposited in Figshare. DOI is listed in the key resource table.

• All original code and commands for all programs have been deposited at Github depository. DOIs are listed in the key resource table. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Sample collection and genomic library preparation

Samples of

METHOD DETAILS

Genome sequencing and assembly

We used a hybrid approach to assemble a draft genome of B. coeca using both longread Oxford Nanopore Technology (ONT) and short-read Illumina sequencing. [START_REF] Miller | Highly Contiguous Genome Assemblies of 15 Drosophila Species Generated Using Nanopore Sequencing[END_REF] Before nanopore sequencing, a size selection was conducted on the DNA using the SRE XS kit from Circulomics (https://www.circulomics.com/, Baltimore, Maryland, USA). The Ligation Sequencing kit SQK-LSK110 from ONT (https://nanoporetech.com/) [START_REF] Lu | Oxford Nanopore MinION Sequencing and Genome Assembly[END_REF] was then used to prepare the samples for nanopore sequencing following the manufacturer's protocol. The library was loaded and sequenced on an R9.4.1 flow cell (ref. FLO-Min106) for sequencing. Raw data were basecalled using Guppy v5.0.11 and the "sup" algorithm.

The ONT raw data size was 4.4 Gb in 1,399,323 reads (mean read length 3,146 kb, longest read of 123.3 Mb), with an N50 of 4,677 kb. Phred scores ranged from 8 to 18, with a median of 13, as assessed by PycoQC.

Estimation of genome size and endosymbionts detection

K-mers frequencies within short-read data were obtained with KMC 3. [START_REF] Kokot | KMC 3: counting and manipulating k-mer statistics[END_REF] Genome size and ploidy were inferred using GenomeScope v2.0 with k-mer size = 21 and Smudgeplot. [START_REF] Ranallo-Benavidez | GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes[END_REF] Contig taxonomy was performed using Blobtools [START_REF] Laetsch | BlobTools: Interrogation of genome assemblies[END_REF] with Diamond as search engine [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF] against the UniProt database using a local copy of the NCBI TaxID file for the taxonomic assignation of the best hit. Minimap2 86 was used for read mapping (Figure S1).

Genome annotation

The B. coeca genome was annotated using Maker v2.31. were obtained from NCBI and used to guide the annotation. Protein-Protein BLAST 2.9.0+ 88 (-evalue 1e-6 -max_hsps 1 -max_target_seqs 1) was then used to assess putative protein functions in B. coeca by comparing the protein sequences given by Maker to the protein sequences from the annotated genome of D. melanogaster. The completeness of genome annotation was assessed using BUSCO at each round and the round with the highest score was retained.

Phylogenomic analysis of the Ephydroidea

Besides our B. coeca assembly, we obtained from NCBI repository genome assemblies for 12 species, transcriptome shotgun assemblies (TSA) for four species, and sequence read runs (SRR) for three species (Table S1). Paired-end DNA raw data of two species, namely Rhinoleucophenga cf. bivisualis and Cacoxenus indagator were assembled using MaSuRCA with default parameters. The transcriptome of Acletoxenus sp. was assembled using Trinity software package 89 on the Galaxy Europe website 90 following standard protocol 91 . BUSCO v.5.0 was used to assess the completeness of those assemblies and to extract single-copy BUSCO genes for all species. Protein sequences of 3,100 single and complete BUSCO genes were aligned using MAFFT [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] and concatenated into a single supermatrix (2,557,349 amino acids). A maximum-likelihood (ML) phylogeny was then inferred for the supermatrix using IqTREE 2 [START_REF] Minh | IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era[END_REF] with 1,000 ultrafast bootstrap iterations [START_REF] Hoang | UFBoot2: improving the ultrafast bootstrap approximation[END_REF] and using the JTT+R substitution model inferred by ModelFinder 95 implemented by IqTree.

Reconstruction of the ancestral ecological niches

For each of the 20 ephydroid species we obtained a predominant ecological niche from the taxonomic literature. [START_REF] Throckmorton | The phylogeny, ecology, and geography of Drosophila[END_REF][START_REF] Markow | Evolutionary genetics of reproductive behavior in Drosophila: connecting the dots[END_REF][START_REF] Ashburner | Entomophagous and other bizarre Drosophilidae[END_REF][START_REF] Okada | Bleeding sap preference of the Drosophilid flies[END_REF] Eight predominant niches were coded as discrete traits, and the Multistate program of the BayesTraits v.4 package 97 was used under the Reverse Jump MCMC model with 1,010,000 chain iterations and a burnin sample of 10,000.

Phylogenomic analysis of Diptera and Hymenoptera

The second phylogenomic analysis involved, besides B. coeca, 25 hymenopteran and 16 dipteran species for which an assembly can be downloaded from the NCBI Genome repository (Table S1). Protein sequences for all species but three, namely Leucophenga varia, Phortica variegata, and Ephydra gracilis, were obtained from NCBI. For these three species, we used the same four-round annotation procedure that we used for B. coeca to identify protein-coding genes and translate their sequences. We used BUSCO to assess the completeness of all annotated and downloaded genomes and their corresponding assemblies. OrthoFinder 46 was used to generate protein sequences of protein-codinggenes of the 42 species and to cluster these sequences into orthogroups. Only the longest isoform (i.e. the primary transcript) was used for genes with multiple isoforms. 79 orthogroups contained a single copy ortholog from each species, and their protein sequences were aligned using MAFFT and concatenated into a single supermatrix (63,192 amino acids). A maximum-likelihood (ML) phylogeny was then inferred for the supermatrix using IqTREE 2 [START_REF] Minh | IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era[END_REF] with 1,000 ultrafast bootstrap iterations [START_REF] Hoang | UFBoot2: improving the ultrafast bootstrap approximation[END_REF] and the JTT+R substitution model inferred by ModelFinder 95 implemented by IqTree.

MCMCTree [START_REF] Yang | PAML 4: phylogenetic analysis by maximum likelihood[END_REF] was used to date the inferred ML trees based on recently published fossil-calibrated phylogenies. First, two time points were obtained for the 42-species phylogeny. These included the divergence between ants and bees between 90-120 myr ago [START_REF] Misof | Phylogenomics resolves the timing and pattern of insect evolution[END_REF] and between Scaptodrosophila and Drosophila between 50-56 myr ago, 100 with a maximum root age for the ancestor of Hymenoptera and Diptera at 344 myr ago. 99

Genome size and gene content evolution

Genome size and gene content (number of OrthoFinder generated protein-codinggenes after retaining the longest isoform for genes with multiple transcripts) inferred for each of the 42 dipteran and hymenopteran genomes were mapped on the phylogenetic tree, and values at the ancestral nodes were inferred and visualized using the fastAnc command in the R package Phytools v0.2.2. 101

Transposons annotation and detection of Horizontal Transposon Transfer (HTT)

Transposons were identified in the 42 dipteran and hymenopteran genomes following a two-step protocol. First, we used RepeatModeler v2.0.1 [START_REF] Flynn | RepeatModeler2 for automated genomic discovery of transposable element families[END_REF] with default parameters to generate a de novo library of repetitive regions. RepeatMasker v 4.0.9 [START_REF] Flynn | RepeatModeler2 for automated genomic discovery of transposable element families[END_REF] was then run with the newly generated library and the options -a (create a .align output file)

and -s (slow search; more sensitive) to create a summary of the families of transposable elements found in each genome along with the percentage of the genome they represent.

To detect possible HTT between Braula coeca and Apis mellifera, we used the B. coeca whole genome as query to perform a blastn similarity search against the whole A.

mellifera genome (all default options, including "-task megablast"). All B. coeca genome regions longer than 299 bp and aligning to A. mellifera with an e-value lower than 0.0001 were extracted and clustered at 80% nucleotide identity threshold with vsearch. 102 The consensus sequence of each of the 50 resulting clusters were used as queries to perform blastx searches on the non-redundant protein database of NCBI using Diamond. [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF] A total of eight consensus sequences had best hits to the Famar1 element previously described in the earwig Forficula auricularia, known to be also present in A. mellifera as a result of horizontal transfer. [START_REF] Barry | A Bacterial Genetic Screen Identifies Functional Coding Sequences of the Insect mariner Transposable Element Famar1 Amplified From the Genome of the Earwig, Forficula auricularia[END_REF]103 To verify that the Famar1-like element from B. coeca has indeed been involved in HTT, we compared the Famar1-like synonymous distance (dS) to a distribution of dS expected under vertical transmission since the last common ancestor of B. coeca and A. mellifera following the approach developed in Zhang et al. 104 This approach assumes that in case of HTT, TE dS should be much lower than dS expected under vertical transmission. Briefly, we calculated the dS over the transposase open reading frame between one copy of the Famar1-like element extracted from C. coeca and another copy of this element from the A. mellifera genome. We then compared this distance to the distribution of dS calculated over 2,179 alignments between single copy BUSCO genes that produce best reciprocal hits in blastp similarity searches. 104 We found that the Famar1-like dS (=0.12) fall below the 0.5% quantile (=1.76) of the distribution of dS calculated for orthologous genes (Figure S2), confirming that the element has been acquired through HTT in B. coeca and A. mellifera. To assess whether the tight ecological interactions existing between B. coeca and A. mellifera might have favored direct transfer of this element between the two species, we assessed how closely related are B. coeca Famar1-like copies to those from A. mellifera. We first screened for the presence of this element in other animal genomes. We used the Famar1 sequence 42 as query to perform online blastn similarity searches (all default options, including "-task megablast") on a total of 8,180 animal genomes belonging to 11 insect orders as well as to Annelida, Chelicerata, Chiroptera, Cnidaria, Myriapoda, Nematoda, Platyhelminthes and Teleostei.

We found full length copies showing >79% of nucleotide identity to this element in a total of 37 species. We aligned up to ten copies from each genome the most similar to Famar1 using Muscle. 105 We then reconstructed a maximum-likelihood phylogeny of these copies using IqTree after nucleotide model detection using ModelFinder. Node support was quantified using ultrafast bootstrap as implemented in IqTree.

Gene family evolution

We used CAFE v. 5 [START_REF] Mendes | CAFE 5 models variation in evolutionary rates among gene families[END_REF] to model and infer gene family evolution. We conducted CAFE5 using an error model on the 19,011 orthgroups generated by OrthoFinder and using the time-calibrated phylogenetic tree of the 42 dipteran and hymenopteran species. Consequently, 5,733 protein-coding genes were assigned to an orthogroup and a gene group. Because some orthogroups can have multiple genes with some assigned to different gene groups, each orthogroup was assigned to the largest gene group of its constituent genes. Orthogroups were then clustered according to their assigned functional groups, e.g., the odorant receptors family contained 16 orthgroups (and 60 D.

melanogaster genes). Because some of the genes found in the orthogroups based on their sequence similarity have no functional annotation in Flybase, the total number of D.

melanogaster protein-coding genes to be grouped into gene groups was 7,820 genes (and 6,317 protein-coding genes for B. coeca).

CAFE5 was then run on the gene groups' gene counts using four different birth rate models in an increasing order (lambda = 1, 2, 3 and 4) and the error model to correct for possible assembly and annotation errors. For each model we run four iterations. The likelihood of only the two simplest models, i.e. one-and two-lambda models, converged across the four iterations. Likelihood ratio test using the lr.test function of the extRemes R package 106 showed that the two-lambda better fit our data. This model imposed a different rate for only B. coeca compared to the rest of the tree and it was chosen for fouriterations of subsequent analyses using the estimated error rate. Multiple testing corrections were conducted using the False Discovery Rate (FDR) analysis of the FDRestimation 108 package implemented in R.

Chemosensory superfamilies evolution

To curate B. coeca gustatory receptors (GR) and odorant receptors (OR) genes, we queried D. melanogaster GRs and ORs protein sequences on B. coeca, L. varia and P.

variegata assemblies using Exonerate ver. 2.2 107 with option -maxintron 2000 and -p pam250. The output, along with the assembly, were fed to InsectOR 54 with option 7tm_7 and 7tm_6 activated for GR and OR analyses, respectively. From the output files, we extracted 300-500 amino acids-long complete sequences with 7tm_7 or 7tm_6 motif detected and with start codon present and no internal stop codon, i.e. pseudogenes excluded. Protein sequences were aligned using MAFFT and a maximum-likelihood phylogenetic tree for each family using IqTREE 2 with the same options as the phylogenomic analysis. The literature was reviewed to classify GRs into bitter, sweet, and CO2 categories 109 and identify volatile ligands eliciting the strongest response in odorant neurons in D. melanogaster. [START_REF] Münch | DoOR 2.0 -Comprehensive mapping of Drosophila melanogaster odorant responses[END_REF] Because CAFE5 inferred ancestral counts for orthogroups with significant deviation only, we estimated and visualized ancestral counts for each orthogroups of these two families using FastAnc command in the R package Phytools. [1]

QUANTIFICATION AND STATISTICAL ANALYSIS

Genome size and protein-coding gene content analyses

D E

[1]

[4]

[0]

[1]

[1]

[1] S3. Gene groups were defined according to D. melanogaster genes clustered with orthologous sequences from 42 dipteran and hymenopteran genomes by OrthoFinder. Puta@ve func@ons of each group are given following FlyBase defini@ons and references therein. Evolu@onary rate was es@mated by CAFE5, with p-values corrected for mul@ple tes@ng using False Discovery Rate (FDR) analysis.

[1] [2] [1] [0] [1] [1] [2] [2]

3 x 10 - 4 ;

 104 FigureS2; TableS1).

  pointing to a possible conserved role in reproduction. Therefore, further functional analyses of the response of candidate ORs to various QMP compounds are required in both D. melanogaster and B. coeca to understand how modifications of these genes in B. coeca might have facilitated the evolution of the bee louse fly inquilinism.

Figure 1 .

 1 Figure 1. The bee louse fly (Braula coeca) is an inquiline of the Western honey bee (Apis mellifera) and has likely evolved from a sap-breeding drosophilid associated with scale insects. See also Table S1. (A) Tens of B. coeca adults preferentially attached to the honey bee queen (© Etienne Minaud). Scale bar = 5 mm. (B) Dorsal view of an adult showing the loss of the wings, halters, and scutum, mesonotum reduction and the legs' robustness. Scale bar = 0.5 mm. (C) Frontal view of an adult showing the reduction of the eyes and the loss of the ocelli. Scale bar = 0.5 mm. (D) Maximum likelihood phylogeny inferred from 3,100 conserved single-copy proteins (2,557,349 amino acids) showing the position of B. coeca (orange) in the subfamily Steganinae (light green) of the Drosophilidae (red). Outgroup species belong to the superfamily Ephydroidea (light blue). All internal nodes had an ultra-fast bootstrap value of 100% except * = 73%. Pie charts at internal nodes indicate the likelihood of ancestral breeding niches inferred from the predominant niches of terminal taxa.

Figure 2 .

 2 Figure 2. Evolution of the bee louse fly inquilinism, its genome size, gene content, and transposable elements in the bee louse fly with evidence for horizontal transfer between the inquiline and its host. See also Table S1 and Figure S2. (A) Fossil-calibrated maximum-likelihood phylogeny inferred from 79 conserved singlecopy proteins (63,192 amino acids) demonstrating major stages in the evolution of the inquiline and its social host. All internal nodes had an ultra-fast bootstrap value of 100 (except when given), with a blue interval indicating a 95% confidence level of divergence time estimate inferred by MCMCTree. The red bar indicates the likely interval of the origin of the bee louse fly-Apis association. Labels 1-6 refer to the major stages mentioned in the text. (B) Genome-size evolution. Red asterisk indicates the estimate for B. coeca. (C) Gene content evolution. Red asterisk indicates the estimate for B. coeca.

Figure 3 .

 3 Figure 3. Evolution of chemosensory receptors gene families in Braula coeca and Drosophila melanogaster. See also Figure S4. (A) Maximum-likelihood phylogeny of gustatory receptors (GRs) with main taste categories color code given in a frame. (B) Maximum-likelihood phylogeny of odorant receptors (ORs) with main ligands for each D. melanogaster receptor given in dark red. L = larva and A = adult expression. For A and B, ultra-fast bootstrap values are given above nodes. Branches are colored according to orthogroups defined by OrthoFinder for 42 dipteran and hymenopteran species. Numbers in broken brackets before each orthogroup reflect the presumed ancestral gene content inferred by phytools.

  Braula coeca were collected from honey bee colonies on the Island of Ouessant in France and kindly provided to us by the Association Conservatoire de l'Abeille Noire Bretonne (A.C.A.N.B.). Genomic DNA was extracted from 15 unsexed individuals conserved in alcohol using the Nucleobond AXG20 kit and buffer set IV from Macherey-Nagel (ref. 740544 and 740604, https://www.mn-net.com, Düren, Germany).

  Figure 2

  -(-)-fenchone (A) propyl acetate (L/A) ethyl trans 2-buteonate (A) 2-methylphenol (L)

Figure 3

 3 Figure 3

  [START_REF] Leger | pycoQC, interactive quality control for Oxford Nanopore Sequencing[END_REF] Illumina paired-end sequencing was performed by Novogene Company Limited (https://en.novogene.com, Cambridge, UK) on the same

	DNA sample. The Illumina sequencing produced 119,719,537 paired 150 bp reads. Phred scores averaged 36 per read as analyzed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used MaSuRCA v4.0.3 diptera_odb10 dataset (C:95.8%[S:94.5%,D:1.3%],F:0.7%,M:3.5%,n:3285 ), and to 93.6% using Merqury.

[START_REF] Zimin | Hybrid assembly of the large and highly repetitive genome of Aegilops tauschii, a progenitor of bread wheat, with the MaSuRCA mega-reads algorithm[END_REF]

to produce the hybrid assembly of our genome using the Cabog assembler. We obtained a final assembly size of 309,35Mb in 2477 contigs, with a N50 of 347,227 bp. The completeness of the assembly was estimated to 95.8% with BUSCO v5.0.0 on the

  10, 33 following Muller et al.'s 87 protocol, wherein multiple rounds of Maker supported by the training of the SNAP v.2006-07-28 34 and Augustus v.3.3.3 35 gene finding and prediction tools, were conducted. RepeatModeler v2.0.1 was first used to identify the repeat-enriched regions that were masked by RepeatMasker v4.0.9 as implemented in Maker. Proteomes of five Drosophila

species, namely D. innubila, D. albomicans, D. bipectinata, D. melanogaster, and D. virilis

  genes are assigned to at least one of 10,670 functional gene groups in the FlyBase database, of which some concern RNA genes that, by definition, are not analyzed by OrthoFinder. Because of the hierarchical nature of the functional gene groups annotation in FlyBase as well as to the pleiotropy of certain genes, each D. melanogaster gene was assigned to its largest group, i.e. the group with biggest number of genes.

The analysis showed that B. coeca has gained 439 while losing 1,517 protein-coding genes in agreement with the low gene content of this species. To gain further functional insights on the ecological or biological relevance of the evolving genes, we grouped OrthoFinder orthogroups into functional gene groups using the customized perl script OG2GG.pl (https://github.com/AmirYassinLab/OG2GG) leveraging the proximity of B. coeca to D. melanogaster. The script assigns each D. melanogaster ortholog to its largest functional gene groups in Flybase 48 ("gene_group_data_fb_2023_02.tsv") and then assigns each orthogroup to the largest gene group of its constituent genes. D. melanogaster has 13,545 protein-coding genes that were clustered into 10,497 orthogroups. However, 8,202 D. melanogaster

Table 1 -Rapidly evolving gene groups in Braula coeca show parallel reduc7on with bees. See also Figure

 1

Acknowledgments

We are very grateful to Marcus Stensmyr for insightful comments on an early draft of the manuscript, the Association Conservatoire de l'Abeille Noire Bretonne (A.C.A.N.B.) for help collecting B. coeca flies, Etienne Minaud for photographic images of B. coeca and Héloïse Muller for help with the genome annotation pipeline, Christian Cheminade and Michael Lang for help translating early literature on Braula. Braula genome sequencing was funded by a grant from Université Paris-Saclay (ADAPAR) to HB. B. pyrosoma B. affinis E. mexicana C. calcarata H. laboriosa O. bicornis M. rotundata N. melanderi M. genalis D. novaeangliae H. anthracinus C. gigas M. pharaonis V. vulgaris N. vitripennis L. boulardi F. arisanus D. alloeum

Report

The genome of the blind bee louse fly reveals deep convergences with its social host and illuminates Drosophila origins Héloïse Bastide 1,* , Hélène Legout 1 , Noé Dogbo 1 , David Ogereau 1 , Carolina Prediger 1 , Julie Carcaud 1 , Jonathan Filée 1 , Lionel Garnery 1 , Clément Gilbert 1 , Frédéric Marion-Poll 1,2 , Fabrice Requier 1 , Jean-Christophe Sandoz 1 , Amir Yassin

Declaration of interests

The authors declare no competing interests.

Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work the authors used Grammarly (Grammarly Inc.) in order to improve language and readability. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

We compared the estimates of genome size and protein-coding gene content of B. coeca to the distributions of these values for the 10 drosophilid genomes (data given in Table S1) using one-sample Student's t test and after testing for normality using Shapiro-Wilk test as implemented in R.

Likelihood Ratio Test (LTR) comparison of CAFE5 models

We compared the likelihood of the two simplest CAFE5 models, i.e. one-and two-lambda models which were the only ones to converge across the four iterations, using the lr.test function of the extRemes R package. 106 For each model, the likelihood estimates were averaged across the four iterations.

False Discovery Rate (FDR) estimation of the p-values of the best CAFE5 model

The two-lambda model had the best likelihood and was therefore subsequently run for four iterations. p-values inferred for each gene group in the branch leading to B. coeca in the iteration with the best likelihood were extracted, and multiple testing corrections were conducted using the False Discovery Rate (FDR) analysis of the FDRestimation 108 package implemented in R.

List of supplementary materials


François Gelis

Grégoire Misguich

Keywords: Out[4]= {KernelObject[1, local, <defunct>], KernelObject[ ...) * .. Parameters .. INTEGER M, N, NRHS PARAMETER Intrinsic Functions .. INTRINSIC INT, MIN

Evolution of 16 confs. on a 32 x 32 x L z grid [2x4 cores 2.93GHz 2x8MB cache] conf. size 1 conf. / 1 core 1 conf. / 4 cores 1 conf. / 8 cores (naive) 1 conf. / 8 cores (improved)

My Laptop [DELL XPS 13 (2015)]

• CPU: 2.4 to 3.0 GHz, 2 cores • Peak performance: 46 GFLOP/sec (double precision) (this would have placed my laptop on the top500 list of the world's most powerful computers, circa 1999...)

Let's do the math...

• At 3.0 GHz, 46 GFLOP/sec is equivalent to 15.3 FLOP/cycle i.e., 7.7 FLOP/cycle/core =⇒ there is a high amount of instruction-level parallelism In this case: each core uses AVX instructions to perform 4 double precision OPs at once, and has two Floating Point Units My Laptop [cont..]

• Memory Bandwidth:

• from RAM: 17-20 GB/sec • from L2 cache: 100-140 GB/sec • from L1 cache: 570-980 GB/sec

• Memory Latency:

• from RAM to CPU: 300 cycles • from L2 cache to CPU: 12 cycles • from L1 cache to CPU: 4 cycles • from L2 (core 1) to L2 (core 2): 90 cycles ←relevant for OpenMP

• When the access pattern is regular, the CPU does some data prefetching, and these latencies are somewhat avoided

• With many cores, the needs for memory bandwidth become more severe

Cache: use correct loop order in nested loops

Vectorization

• The compiler tries automatically with -O3, but not perfect

• Biggest obstacle: dependence among loop indices:

for (int j = 1; j < N; ++j){ result[j] = 1.0/(1.0+result[j-1]); }
• Less if the dependence has a range longer than the vector length Parallelization on shared memory:

for (int j = 4; j < N; ++j){ result[j] = 1.0/(1.0+result[j- 4
• Very easy to implement via OpenMP, with tiny code modifications

• Should scale perfectly for independent tasks Large scale eigenvalue problems. For (usually sparse) nxn matrices with n which can be as large as 10^8, or even more ! Ex: iterative algo. to find the largest eigenvalue, without storing the matrix M (just provide v→Mv).

Can be used from Python/SciPy Examples for several LAPACK routines that solve systems of linear equations.

Linear Least Squares Problems

Examples for some of the LAPACK routines that find solutions to linear least squares problems.

Symmetric Eigenproblems

Symmetric Eigenproblems has examples for LAPACK routines that compute eigenvalues and eigenvectors of real symmetric and complex Hermitian matrices.

Nonsymmetric Eigenproblems

Nonsymmetic Eigenproblems provides examples for ?geev, one of several LAPACK routines that compute eigenvalues and eigenvectors of general matrices.

Singular Value Decomposition

Examples for LAPACK routines that compute the singular value decomposition of a general rectangular matrix. Compilation and output :

Lapack-MKL example in Fortran

First call to DGELSD with LWORK=-1 → LAPACK returns the optimal size LWORK of the workspace array WORK.

Actual calculation

Lapack-MKL example in Fortran checking what linear algebra library that is actually used $ ldd diag.exe linux-vdso.so.1 => (0x00007ffeddd7c000) libarmadillo.so.9 => /usr/local/install/armadillo-9.200.6/lib64/libarmadillo.so.9 (0x00007f2c30549000) libstdc++.so.6 => /usr/local/install/gcc-4.8.0/lib64/libstdc++.so.6 (0x00007f2c3023f000) libm.so.6 => /lib64/libm.so.6 (0x000000338f200000) libgcc_s.so.1 => /usr/local/install/gcc-4.8.0/lib64/libgcc_s.so.1 (0x00007f2c30010000) libc.so.6 => /lib64/libc.so.6 (0x000000338e200000) libmkl_rt.so => /opt/intel/composer_xe_2013.5.192/mkl/lib/intel64/libmkl_rt.so (0x00007f2c2fb02000) libhdf5.so.6 => /usr/lib64/libhdf5.so.6 (0x0000003390a00000) libz.so.1 => /lib64/libz.so.1 (0x000000338ee00000) /lib64/ld-linux-x86-64.so.2 (0x000000338de00000) libdl.so.2 => /lib64/libdl.so.2 (0x000000338ea00000)

MKL used here The corresponding numpy will be using the MKL lib.

Process

• Each process has its own memory space. → they are somewhat independent • Switching from one process requires some interaction with the operating system → slow switching Thread

• Threads (associated to a given process) share the same memory space → Threads can share information easily /quickly • There is no memory "protection" between the threads of the same process → responsibility of the programmer • Threads have little information of their own → faster to create than processes Can be useful for some I/O tasks (because the CPU will be waiting for some remote server, etc.), but not really for computations. … will not be discussed here.

• Multiprocessing library (import multiprocessing) Allows to perform tasks simultaneously (using processes instead of threads)

We will present a few examples using: For loops Schedule(static):

• Process, Queue • Pool, map, imap
1 2 3 … thread #0 thread #1 thread #2 thread #3 thread #4 thread #5 33
Example above: #pragma omp parallel num_threads(6) { #pragma omp for schedule(static) for (int i=0;i<33;i++) … }

• Iterations are divided into 'chunks' of size chunk_size and distributed cyclically to the threads. • If the chunk_size is not specified, the iterations are divided into (almost) equal chunks, and each thread executes one chunk (example below).

For loops Schedule(static)

Another example: #pragma omp parallel num_threads(6) { #pragma omp for schedule(static,2) for (int i=0;i<33;i++) … }

1 2 3 … 33 thread #0 thread #1 thread #2 thread #3 thread #4 thread #5 thread #0 thread #1 thread #2 thread #3 thread #4 thread #5 thread #0 thread #1 thread #2 thread #3 thread #4
Advantage of large chunks: less overhead, cache friendly Advantage of small chunks: better load balance between the threads

For loops

Schedule(dynamic)

• The iterations are divided in chunks of size chunk_size (the last one can be smaller). • When a thread is idle, it is assigned a new chunk (first come, first served).

• If chunk_size is not specified, it is set to 1. 1 2 3 … thread #0 thread #1
1 and 3: long tasks.

For loops Schedule(guided)

• Similar to dynamic, but the chunk size is initially large, and decreased gradually • The size of a chunk is proportional to the number of remaining iterations, divide by the number of threads) • chunk_size specifies the minimum size of the chunks. If not specified, this minimal size is set to 1.

For loops

Schedule(runtime)

[0]=0 A[1]=0 A[2]=0 A[3]=0 A[4]=0 A[5]=1 A[6]=1 A[7]=1 A[8]=0 A[9]=1

Message Passing Interface (MPI)

• MPI provides high level functions to exchange data between jobs on several nodes, that hide the network details

Variants of MPI_Send

• MPI_Send : blocks until the array containing the data to be sent out can be reused safely (does not imply that the data has already reached its destination, since it could be buffered)

• This function has several variants

Variants of MPI_Recv

• MPI_Recv : blocks until the array in which the data to be received arrives is ready to be used

• This function has one variant Easiest: independent computations on each node

• Each MPI job performs a (lengthy) computation, independently of the other jobs

• These computations are parallelized using OpenMP

• Output is sent back to MASTER node for final processing (e.g., averaging, saving in a file, etc...). For this, the MASTER node forks a special thread that collects this output • Input data:

• Monte-Carlo: input is a random configuration, generated locally on each node (make sure each MPI job uses its own RNG seed)

• Alternative: the MASTER node sends tasks taken from a list to each node as soon as a result returns

Threads on MASTER node

• Receiving/processing computed results is a very light task compared to the computations ⇒ the MASTER node can also perform computations without significant penalty

• Call omp_set_nested(1) on all nodes

• Create two parallel sections:

• First section: executed only by MASTER (empty on other nodes) Post MPI_Recv to receive results, post-process results, MPI_Send to send tasks, etc...

• Second section: executed by everybody performs actual computations ⇒ further fork in Ncores threads

Harder: Two-level slicing of the computation domain

• Large array to be processed (e.g., evolved in time)

• Divide the array in N cores × N threads slices, and assign a slice to each thread

• Case 1: evolution is "almost" local, i.e. depends at most on a few neighboring sites (e.g., discretization of a Laplacian)

At the beginning of each timestep, each thread must be given a copy of the layers just before and just after its slice

• Case 2: evolution is completely non-local (update of a point i depends on all other points) 2.a. Array is small enough: each node can have its own copy of the full array (it must be refreshed at the beginning of each timestep) 2.b. Array is too large: communications will probably make parallelization very inefficient

Boltzmann Equation

• Two-body elastic interactions

• Spatially homogenous

• Isotropic particle distribution

• Scattering amplitude may be momentum dependent

∂ t f 1 = 1 E 1 p 2,3,4
M(1, 2, 3, 4) 2 f 3 f 4 (1+f 1 )(1+f 2 )-f 1 f 2 (1+f 3 )(1+f 4 )

• Collision integral reduces to 4-dim integral thanks to momentum conservation and isotropy

3. 0

 0 GHz × (2 cores) × (2 FPUs) × 4(AVX) = 48 GFLOP/sec • What bandwidth do we need for this? To simplify, assume 1 FLOP = 2 reads + 1 write 48 GFLOP/sec =48 × 8 (Bytes in a double) × 3 (2 reads + 1 write) GB/sec =1152 GB/sec François Gelis, Grégoire Misguich, June 2019

  for (i=0;i<N;i++){ for (j=0;j<N;j++){ b[j+N*i] = a[i+N*j]; } } • Improvement: slice the j loop in blocks of size B: for (jb=0;jb<N;jb+=B){ for (i=0;i<N;i++){ for (j=0;j<B;j++){ b[j+jb+N*i] = a[i+N*(j+jb)]; } } } Cache: sometimes, "loop tiling" helps • How does it work? • This transformation makes the two innermost loops work on a contiguous range of size N * B doubles • For N = 2 13 , the optimal block-size is around B = 2 4 (empirical) • 2 13 × 2 4 × 8 (Bytes in a double) = 1 MByte = order of L2 cache• Note: the compiler can do this transformation automatically (but perhaps not choose the best block-size)

×

  (# FPUs) × (Width of vector inst.) well exploited only for appropriate code × (# cores) try this?

  's documentation] « Maple provides tools for two different types of parallel programming. The Task Programming Model enables parallelism by executing multiple tasks within a single process. The second type of parallelism comes from the Grid package, which enables parallelism by starting multiple processes."Remark: the task model is somewhat analogous Threads/OpenMP, and Grid is analogous to MPI. Both OpenMP and MPI will be presented in these lectures2. Automatic parallelizationfinish Shared Memory & Synchronization » SetSharedVariable -specify symbols with values to synchronize across subkernels SetSharedFunction -specify functions whose evaluations are to be synchronized $SharedVariables ▪ $SharedFunctions ▪ UnsetShared ▪ CriticalSectio n Setup and Configuration » LaunchKernels -launch a specified number of subkernels $KernelCount -number of running subkernels $KernelID ▪ Kernels ▪ AbortKernels ▪ CloseKernels ▪ ... $ProcessorCount -number of processor cores on the current computer Multi-Processor and Multicore Computation Compile -create compiled functions that run in parallel Parallelization -execute compiled functions in parallel CompilationTarget -create machine-level parallel compiled functions GPU Computing » CUDAFunctionLoad -load a function to run on a GPU using CUDA OpenCLFunctionLoad -load a function to run on a GPU using OpenCL File-Based Parallelism FileSystemScan ▪ FileSystemMap 1. Introduction & hardware aspects (FG) 2. A few words about Maple & Mathematica 3Here: a few simple examples showing how to call some parallel linear algebra libraries in numerical calculations (numerical) Linear algebra • Basic Linear Algebra Subroutines: BLAS • vector op. (=level 1) • matrix-vector (=level 2) • matrix-matrix mult. & triangular inversion (=level 3) • Many implementations but standardized interface • Discussed here: Intel MKL & OpenBlas (multi-threaded = parallelized for shared-memory architectures) • More advanced operations Linear Algebra Package: LAPACK ('90, Fortran 77) • Call the BLAS routines • Matrix diagonalization, linear systems and eigenvalue problems • Matrix decompositions: LU, QR, SVD, Cholesky • Many implementations Used in most scientific softwares & libraries (Python/Numpy, Maple, Mathematica, Matlab, …) A few other useful libs … for BIG matrices • ARPACK =Implicitly Restarted Arnoldi Method (~Lanczos for Hermitian cases)

  for(j = 0; j<N; j++ ) a[i*LDA+j].real=rand(),a[i*LDA+j].imag=rand( the num. of threads Call to Lapack. This version takes care of the workspace memory management (contrary to the Fortran version) w: array containing the eigen. vals. for(j = 0; j<N; j++ ) a[i*LDA+j]=lapack_make_complex_double(rand() ,rand()); // LAPACKE_zheevd: computes all // eigenvalues and eigenvectors of a // complex Hermitian matrix A using divide // and conquer algorithm info = LAPACKE_zheevd( LAPACK_ROW_MAJOR, 'V', 'L', N, a, LDA, w ); /* Check for convergence */ if( info > 0 ) { printf( "The algorithm failed to compute eigenvalues.\n" ); exit( 1 ); } /* Print the extreme eigenvalues */ printf("Smallest eigen value=%6.2f\n",w[0]); printf("Largest value=%6.2f\n\n",w[N-1]); } OpenBLAS/CBLAS dgemm example in C #include <stdio.h> #include <stdlib.h> #include <cblas.h> int main() { int N=10000,N2,i,j; N2=N*N; //Memory allocation for the arrays: double *A, *B, *C; A = (double *)malloc( N2*sizeof( double ) ); B = (double *)malloc( N2*sizeof( double ) ); C = (double *)malloc( N2*sizeof( double ) ); for (i = 0; i < (N2); i++) A[i] = (double)(i+1), B[i] = (double)(-i-1), C[i] = 0.0; printf ("Computing matrix product using OpenBLAS dgemm function via CBLAS interface ...\n"); cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N, N, N, 1.0, A, N, B, N, 1.0, C, N); printf ("done.\n\n"); return 0; } Matrix-matrix multiplcation (BLAS) of double precision general matrices Specify the « order » of the matrix elements in memory: A[i][j]=A[j+LDA*i] row-major (C/C++ order) A[i][j]=A[i+LDA*j] column-major (=Fortran order) Compilation of the OpenBLAS examples -use of make & makefile #specify below the path to the OpenBLAS library files (like libopenblas.a and cblas.h) OPEN_BLAS_LIB= ../OpenBLAS #specify below the path to the lapacke.h header file LAPACKE_INC= ../OpenBLAS/lapack-netlib/LAPACKE/include dgemm_example_mini.exe: dgemm_example_mini.c gcc $< -o $@ -L $(OPEN_BLAS_LIB)-I $(OPEN_BLAS_LIB) -lopenblas -pthread zheevd_example.exe: zheevd_example.c gcc $< -o $@ -I $(LAPACKE_INC) -L $(OPEN_BLAS_LIB) -lopenblas -fopenmp -lgfortran clean: \rm *.exe all: dgemm_example.exe zheevd_example.exe To compile: make dgemm_example.exe or make all or make -j 2 all -j option: use several threads/cores to compile multiple files in parallel makefile $< : 1st pre-requiste (usually the source file) $@ : target (executable name) $ make dgemm_example.exe gcc dgemm_example.c -o dgemm_example.exe -L ../OpenBLAS -I ../OpenBLAS -lopenblas -fopenmp -lrt $ export OMP_NUM_THREADS=1; ./dgemm_example.exe 2000 Initializing the matrices ... done. Computing matrix product using OpenBLAS dgemm function via CBLAS interface... done. Elaspe time (s): 0.719543 22.2308 GFlops $ export OMP_NUM_THREADS=10; ./dgemm_example.exe 2000 Initializing the matrices ... done. Computing matrix product using OpenBLAS dgemm function via CBLAS interface... done. Elaspe time (s): 0.0814542 196.38 GFlops Close to the peak power of the CPU (here Xeon E5-2630 v2 @ 2.60GHz / 15360 KB Cache). Check with cat /proc/cpuinfo ) compilation ARMA_INC=/usr/local/install/armadillo-9.200.6/include ARMA_LIB=/usr/local/install/armadillo-9.200.6/lib64/ diag.exe: diag.cpp g++ -std=gnu++11 $< -o $@ -I $(ARMA_INC) -L $(ARMA_LIB) -larmadillo mult.exe: mult.cpp g++ -std=gnu++11 $< -o $@ -I $(ARMA_INC) -L $(ARMA_LIB) -larmadillo computing

  BLAS & Lapack from Python/Numpy Python → Numpy → Linalg → LAPACK → BLAS Different possible implementations • standard BLAS/LAPACK (Netlib) • ATLAS (Automatically Tuned Linear Algebra Software) • OpenBLAS • MKL • … Multi-threaded /parallel Python/Numpy linalg Check the version of Lapack & BLAS numpy is linked to: >>> import numpy as np >>> np.__config__.show() Python 2.7.5 (default, Mar 20 2015, 15:33:03) [GCC 4.4.7 20120313 (Red Hat 4.4.7-11)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import numpy as np >>> np.__config__.show() lapack_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('HAVE_CBLAS', None)] language = c blas_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('HAVE_CBLAS', None)] language = c openblas_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('HAVE_CBLAS', None)] language = c (. . .) Python/Numpy linalg What if my numpy version is not linked to a parallel linear algebra lib. ? • Install OpenBLAS: > sudo apt-get install libopenblas-dev (this will hopefully replace the previous BLAS lib. by OpenBLAS in Numpy) • Or the Intel Python distribution https://software.intel.com/en-us/articles/installingintel-free-libs-and-python-apt-repo (…) > sudo apt-get install intelpython3

•#

  Matrix diagonalization exampleSpecify the #of threads using an environnent variable (bash shell): > export OMP_NUM_THREADS=4 Python/Numpy linalg import numpy as np import numpy.random as npr import time npr.seed(2019) n=3000 A = npr.randn(n,n) t = time.time() v = np.linalg.eigvals(A) td = time.time()t print(" Time=%0.4f s" % (td)) 1. Introduction & hardware aspects (FG) 2. A few words about Maple & Mathematica 3FFTW: a high performance implementation (in C) Fastest Fourier transform in the West • http://www.fftw.org/ • Multi-threaded • portable • Open source (GPL license) • Can be used from Python (pyFFTW) (discrete) Fourier transform example using multi-threaded FFTW v3 A given application/program may run different processes and different threads

•

  Thread library (import threading) Can start several threads, but they will not run simultaneously (because of the Global Interpreter Lock -a.k.a. GIL).

•'

  import os, math from multiprocessing import Process,Queue def my_function(r,q): proc = os.getpid() i1=r[0]**3 i2=r[1]**3 print('Process #{0} will sum from {1} to {2}'.format(proc,i1,i2)) sum=0.0 for i in xrange (i1,i2): sum+=math.sin(i) q.put((i1,i2,sum)) ranges = [ [1,100], [201,300], [300,400],[401,500] ] list_of_procs = [] q=Queue() for r in ranges: p = Process(target=my_function, args=(r,q)) list_of_procs.append(p) p.start() for p in list_of_procs: p.join() results=[q.get() for p in list_of_procs] print(results) q: Object where each process can write its result Processes do not share memory, which means that the global variables are copied, hence their value in the original process do not change. This example does not work (i.e. output is [0,0,0,0]) import os, math from multiprocessing import Process A=[0,0,0,0] def my_function(i): A[i]=A[i]+1 return list_of_procs = [] for r in range(4): p = Process(target=my_function, args=(r,)) list_of_procs.append(p) p.start() for p in list_of_procs: p.join() print(A) Python/Multiprocessing Multiprocessing.Process linear algebra: processes versus threads Single Python processe using ~24 threads (and 24 CPU cores) import multiprocessing as mp from sympy import * p=mp.Pool(processes=3) x = symbols('x') print(p.map(integrate,[x,x**2,x**3])) Python/multiprocessing Multiprocessing.Pool & Sympy import multiprocessing as mp from sympy import * import sys x = symbols('x') A= [ exp(x),sin(x),cos(x),cosh(x),sinh(x)] def my_func(f): count=0 while (count<10000) : f=diff(f,x) count=count+1 return f # Pass the wanted number of process as commandline argument np=int(sys.argv[1]) p=mp.Pool(processes=np) results=p.map(my_func,A) print(results) Based on threads • For shared-memory architectures • Standardized • Mature (goes back to the 90's) • Portable (supported by many compilers, systems and languages) compiler directives to your existing sequential code (written in C, C++ or Fortran) to tell: • Which instructions should executed in parallel • How to distribute (and synchronize) the instructions over the threads • How to distribute/share the data over the threads OpenMP Hello world' with OpenMP Fotran example. Compile with > ifort -fopenmp hello.f90 program hello USE OMP_LIB PRINT *,"Hello, I am thread #", OMP_GET_THREAD_NUM() PRINT *,"a.k.a. the master thread" !$OMP PARALLEL NUM_THREADS(3) PRINT *,"I am thread #", OMP_GET_THREAD_NUM() !$OMP END PARALLEL end program hello OpenMP compiler directives Parallel section =creation of a team of threads [misguich@totoro OpenMP]$ ./a.out Hello, I am thread # 0 a.k.a. the master thread I am thread # 0 Implicit barrier at the end of the parallel section (wait until all the threads in the team have reach this point) 'Hello world' with OpenMP C++ example. Compile with: icc -fopenmp hello.cpp Output: #include <stdio.h> #include <omp.h> int main() { #pragma omp parallel num_threads(3) { printf("Hello, I am thread %d/%d\n",omp_get_thread_num(),omp_get_num_threads()); } } ./hello-cpp.exe Hello, I am thread 0/3 Hello, I am thread 2/3 Hello, I am thread 1/3 For loops with reduction in Fortran program pi USE OMP_LIB INTEGER n, i DOUBLE PRECISION sum, x n=1e9 !$OMP PARALLEL DO REDUCTION(+:sum) private(x) do i=0,n-1 x=(i*1.0)/n sum=sum+sqrt(1-x*x) enddo !$OMP END PARALLEL DO print *, sum/n*4 end program pi $make pi_f90.exe gfortran -fopenmp pi_f90.f90 -o pi_f90.exe $ export OMP_NUM_THREADS=10;time ./pi_f90.exe

3

  Example of task submission script (for PBS/Torque)Listing 1: script.pbs #PBS -S /bin/bash #PBS -N boltzmann #PBS -e job.err #PBS -o job.log #PBS -m abe #PBS -M francois.gelis@ipht.fr #PBS -l nodes=32:ppn=16 module load openmpi/1.6.4 cd $PBS_O_WORKDIR mpirun -npernode 1 ./my_program • This example will start 32 copies of my_program (one per node) • Then, do: qsub script.pbs • Other commands: qstat, qdel • This is sufficient to start independent tasks on several nodes • Non-interactive: I/O to files only François Gelis, Grégoire Misguich, June 2019

  Collective communication: Broadcast• MPI_Bcast(Buf,N,Type,0,MPI_COMM_WORLD) François Gelis, Grégoire Misguich, June 2019 Collective communication: Gather • MPI_Gather(Source,N,Type,Dest,N,Type,0,MPI_COMM_WORLD) • Variant: MPI_Gatherv: gather variable size chunks of data, and place them at variable offsets in the Dest array Collective communication: AllGather • MPI_Allgather(Source,N,Type,Dest,N,Type,MPI_COMM_WORLD) François Gelis, Grégoire Misguich, June 2019 • Superimposing MPI and OpenMP parallelization poses no problem • In general, each MPI job should fork in a number of threads equal to the number of physical cores on one node • Various hybrid scenarios are possible François Gelis, Grégoire Misguich, June 2019

out DGELSD Example Program Results Minimum norm solution -0.69 -0.24 0.06 -0.80 -0.08

	X :	S	REAL
		D	DOUBLE PRECISION
		C	COMPLEX
		Z	COMPLEX*16 or DOUBLE COMPLEX
	YY:	BD DI GB GE	bidiagonal diagonal general band general	0.21 0.12 -0.65 0.29 -0.24 0.38 0.42 0.29 0.35 -0.30
		GG	general matrices, generalized problem (i.e., a pair of general
	matrices)			Effective rank =	4
		GT (…)	general tridiagonal	Singular values 18.66 15.99 10.01	8.51

DGELSD/from Intel's website (part 1/3) Here minimize 𝑏 -𝐴 Ԧ 𝑥 2 with 𝐴 a rectangular matrix and 𝑏 a vector.

Naming convention of the LAPACK routines : XYYZZZ ZZZ: Type of computation. Here LSD stands for minimum norm solution to a linear least squares problem using the singular value decomposition of A and a divide and conquer $ ifort -mkl DGELSD_example.f $ ./a.

•

  The scheduling method is decided only during the execution (=runtime), according to the environment variable OMP_SCHEDULE (or using omp_set_schedule(…));

	#include <stdio.h>
	#include <omp.h>	$ export OMP_SCHEDULE=guided,1 $ ./for-schedule.exe
	int main() {	A
	int n=6;int A[n];
	#pragma omp parallel num_threads(2)
	{
	int id=omp_get_thread_num();
	#pragma omp for schedule(runtime)
	for (int i=0;i<n;i++) A[i]=id;
	}
	for (int i=0;i<n;i++)
	printf("A[%d]=%d\n",i,A[i]);
	}

  LockApplications: computing the number of neighbors of a given node in a network

				2 2
				3
	for (i=0; i<Nv; i++) omp_init_lock(&locks[i]); for (i=0; i<Nv; i++) omp_init_lock(&locks[i]);	BUG !	1 1	0	3 3
	#pragma omp parallel for #pragma omp parallel for		2	4 4
	for (j=0; j<Nb j++){ for (j=0; j<Nb j++){
	omp_set_lock(&locks[bondA[j]]); omp_set_lock(&locks[bondA[j]]);			1
	omp_set_lock(&locks[bondB[j]]); degree[bondA[j]]++;		0 0
	degree[bondA[j]]++; omp_unset_lock(&locks[bondA[j]]); degree[bondB[j]]++; omp_unset_lock(&locks[bondA[j]]); omp_set_lock(&locks[bondB[j]]); degree[bondB[j]]++;		Nv: number of nodes=5 Nv: number of nodes=5 Nb: number of bonds=5 bondA[0]=1;bondsB[0]=4; bondA[1]=4;bondsB[1]=0; Nb: number of bonds=5 bondA[0]=1;bondsB[0]=4; bondA[1]=4;bondsB[1]=0;
	omp_unset_lock(&locks[bondB[j]]); omp_unset_lock(&locks[bondB[j]]);		bondA[2]=0;bondsB[2]=1; bondA[2]=0;bondsB[2]=1;
	} }		bondA[3]=3;bondsB[3]=2; bondA[3]=3;bondsB[3]=2;
			bondA[4]=1;bondsB[4]=1; bondA[4]=1;bondsB[4]=1;

François Gelis, Grégoire Misguich, June 2019

Example of (symmetric) matrix diagonalizationCan call OpenBLAS or MKL #include <iostream> #include <armadillo> using namespace arma; int main() { const int N=5000; size_t dim=N; mat A(dim, dim, arma::fill::randu); vec eigval; mat eigvec; eig_sym(eigval,eigvec,A); cout<<"1st eigenvalue="<<eigval[0]<<"\tLast="<<eigval[dim-1]<<endl; return 0;

Matrix-matrix multiplication#include <iostream> #include <armadillo> using namespace arma; int main() { const int N=10000; size_t dim=N; mat A(dim, dim, arma::fill::randu); mat B(dim, dim, arma::fill::randu); mat C=A*B; return 0; }

workers (=processes), each one using ~12 threads (and 12 CPU cores) [mutli-threaded linear algebra lib.]

Multiprocessing.Pool &

Multiprocessing.Pool & Sympy

my_file.write(str(n)+"\t"+str(td)+"\n") Eigenvals of (10,10) matrix in 0.0008 s Eigenvals of (20,20) matrix in 0.0002 s Eigenvals of (50,50) matrix in 0.0006 s Eigenvals of (100,100) matrix in 0.0359 s Eigenvals of (200,200) matrix in 0.0393 s Eigenvals of (300,300) matrix in 0.0996 s Eigenvals of (400,400) matrix in 0.1522 s Eigenvals of (500,500) matrix in 0.2451 s Eigenvals of (600,600) matrix in 0.3967 s Eigenvals of (700,700) matrix in 0.4587 s Eigenvals of (800,800) matrix in 0.5750 s Eigenvals of (900,900) matrix in 0.6680 s Eigenvals of (1000,1000) matrix in 0.7734 s Eigenvals of (2000,2000) matrix in 2.4721 s Eigenvals of (3000,3000) matrix in 7.4194 s Eigenvals of (4000,4000) matrix in 14.2675 s Eigenvals of (6000,6000) matrix in 39.1304 s Eigenvals of (8000,8000) matrix in 68.1390 s

shared, private

Inside a parallel region:

• Shared variables can be read and written by all the threads. Be careful with potential race conditions. If two threads simultaneously write at the same memory location (variable), or if a threads reads it while another one writes on it, the result is potentially random (possibility of corrupted data). There will be no error message !

• If a variable is private, each thread has its own copy. If a variable existed with the same name before the parallel construct, it is not affected when exiting the parallel region. • By default, variables declared outside the parallel regions are shared, and those declared inside are private. • When entering a parallel region, the private variables are not initialized. In C++ they are created using the default constructor firstprivate, lastprivate

• Firstprivate: special case of private variable, where each local copy is initialized from the value of the variable with the same name before the beginning of the parallel region

• Lastprivate: special case of private variable for parallel section or parallel for, where, at the end of the parallel region, the variable with the same name outside the parallel region gets the value of local copy of the thread doing the last iteration (or last section). printf("task2 done.\n"); } } printf("sum=%g\n",sum); return 0; } Required, to ensure that threads do not attempt to update the shared variable sum simultaneously.

Critical and atomic

Intructions or blocks which must be executed one thread at a time sum is a shared variable before the parallel section. In the parallel for loop a private copy of sum is created for each thread. At the end of the loop the private copies are combined using the operation '+'. I absence of the nowait option, there is an implicit barrier at the end of the loop. All the threads will wait that all the loop iterations are completed before going on. Same implicit barrier at the end of sections or single directives. Here the use of atomic would have not been possible, since vector::push_back(…) is not an "atomic" statement. critical would have been possible, but slower.

Tasks

Useful to parallelized « irregular » problems, unbounded loops, or recursive algorithms. (since OpenMP3).

• Each time a thread reaches a task directive, the corresponding unit of work is added to a queue, and that thread can continue. • A thread of the team (the same or another one) will execute the task (now of later). • All the tasks created by any thread in the current team will be completed before exiting the parallel region. Cutoff needed for performance, to avoid creating too many very small tasks.

Sort the list

wait that both tasks above are finished before returning the result i+j

Only one thread of the team calls fibonacci.

Recursive calls → exploration of a binary tree

Distributed memory "computers"

• Organization of a computer cluster:

• one front-end node for compilation and administration tasks • many nodes for computations, not directly accessible • computation tasks submitted via a batch system • The nodes are connected via a network • Types of network connections:

Hybrid: MPI+OpenMP

François Gelis, Grégoire Misguich, June 2019

Case study: deterministic Boltzmann solver

François Gelis, Grégoire Misguich, June 2019

Sketch of the Algorithm

• At each time-step:

1. For each i, compute ∆f i (given by a 4-dim integral ⇒ slow)

If FALSE, re-run the previous loop with a smaller timestep 3. Do f i + ∆f i → f i and return to step 1

• Note: cost of computing ∆f i not uniform (50% variation) Thus, we expect that the parallelization of the "big loop" will not be perfectly efficient (the computation time will align to that of the slowest bins)

Sequential version (sketch of the relevant bit of code)

// Core of the function that evolves f[ i ]

double * df = (double * )malloc(N * sizeof(double)); for ( i=0;i<N;i++) df[ i ] = dt * C(i , f ) ; // depends on f[k] with k!=i ; computation of C(i,f) very slow


Khalil Aliane
email: khalil.aliane@etu.univ-amu.fr

Damien F G Minenna
email: damien.minenna@protonmail.com

Yves Elskens
email: yves.elskens@univ-amu.fr

Alexandre Frederic Andre
email: frederic.andre@thalesgroup.com

Fr Éd Éric Andr É B

Alexandre Poy
email: alexandre.poye@univ-amu.fr

Time-Domain non-linear simulations of a 3 meter long traveling wave tube

Keywords: metric traveling-wave tubes TWTs, DIMOHA, simulation, time-domain, AM/AM, AM/PM, wave-particle interaction, Hamiltonian dynamics, N-body dynamics, large-signal amplifiers, non linear regime, trapping

come

INTRODUCTION

We presented previously our new time-domain Hamiltonian model for TWTs and similar vacuum structures with a periodic slow wave structure (SWS), which relies on two main features, namely a model order reduction for the electromagnetic wave representations and an N -body Hamiltonian to describe the beam-wave interaction. The order reduction enables our algorithm DIMOHA to be much faster than alternative PIC codes despite the big number involved in an N -body description. One of DIMOHA's advantages is that, since it relies on a self-consistent Hamiltonian, it fully reproduces the physics in the SWS, including non linear phenomena, such as intermodulation, trapping etc. . . Those features were successfully assessed against measurement for several industrial TWTs. However, such devices are primarily designed for broadband or high power amplification [4], and thus do not reach far beyond saturation. One way to assess DIMOHA's ability to investigate strong non linear behaviours is to apply it to meter long, low frequency TWTs. Though uncommon, such TWTs have been used [1,2,3] primarily as investigation tools for beamplasma instabilities, since they share mathematically equivalent evolution equations. Unlike their industrial centimeter long counterparts, meter long TWTs can go deep in the non-linear regime, where several trapping oscillations are observed.

One such TWT is the San Diego university 3 meter long tube, built specificaly to study the beam-plasma instabilty. This TWT is 3 to 4 times longer than commercial TWTs in terms of wavelengths, and contains a continuous probing system, allowing to follow the RF power all along the delay line [1,2]. Another metric TWT is the 4 meter long tube built and developed in the PIIM laboratory at Aix-Marseille university [3], which contains 4 moving probes providing also access to the power along the TWT. Here we compare DIMOHA's simulations of the 3 meter long San Diego TWT, studies of the 4 meter long tube are undergoing.

ALGORITHM

DIMOHA is a 1-D specialized model in time domain and large-signal regime, based on an N -body (a.k.a. many body) Hamiltonian description [5,6,[START_REF] André | Proc. 14th Int. Vac. Electron. Conf. (IVEC)[END_REF][START_REF] André | Proc. 16th Int. Vac. Electron. Conf. (IVEC)[END_REF][START_REF] Minenna | Proc. 18th Int. Vac. Electron. Conf. (IVEC)[END_REF][START_REF] Minenna | Proc. 20th Int. Vac. Electron. Conf. (IVEC)[END_REF]. This symplectic algorithm [START_REF] Minenna | [END_REF] is written in modern Fortran and parallelized with the Message Passing Interface (MPI) libraries. DIMOHA bears several advantages, including a its grand versatility, as it supports any periodic slow wave structure and as it has been tested and validated for a large range of industrial TWTs of different geometries and types [START_REF] André | Proc. 21st Int. Vac. Electron. Conf. (IVEC)[END_REF][START_REF] André | [END_REF], and a relatively small computing cost as well as a small number of input parameters requirement. It is the discrete method used in our model that enables DI-MOHA to be much faster than PIC codes : electromagnetic fields are decomposed on space cells of the size of the helix pitch (more generally, the length of the periodic structure cell). Even for a several meter long TWT, comprising thousands of pitches, we only need to adjust the numerical parameters of our algorithm, and optimise the parallelization and memory management to successfully simulate meter long TWT reasonably fast.

VALIDATION

We use DIMOHA to simulate the electromagnetic power along the delay line in the San Diego 3 meter long TWT at 195 MHz. In figure 1, and we compare our simulation (in blue) to San Diego observations (in orange) [1,2]. We obtain a good agreement in the linear regime. After the saturation point, discrepancies grow bigger : while we observe the same number of trapping oscillations, their measured period is a little longer than what DIMOHA predicts. This can be due to several reasons, the most likely one being an improperly adapted space charge model. DIMOHA also enables us to follow the electron positions and velocities at any given time of the simulation (figure 2).

CONCLUSION

While previously tested and validated for several industrial TWTs, DIMOHA proves to be also a promising tool for more research oriented plasma-beam physics and deep non linear regime investigation. Currently, we are testing our algorithm on the PIIM laboratory homemade 4 meter long TWT.

Figure 1 .

 1 Figure 1. Time average RF power at F = 195 MHz as function of the longitudinal position inside a 3 m long TWT at time t = 0.7µs.

Figure 2 .

 2 Figure 2. Zoom at the macro-electrons positions and velocities (blue dots), normalized by the beam velocity v 0 , inside the 3 m long TWT at time t = 0.7µs simulated with DIMOHA at f = 195 MHz. The dashed line corresponds to the phase velocity v φ /v 0 . Electrons slower than their initial speed have lost their momenta to the benefit of the wave.